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I.  INTRODUCTION 


During  the  contract  period  the  work  focused  on  two  activities:  studies  of  high  energy 
electron  acceleration  and  studies  of  the  Alfven  maser.  The  work  on  acceleration  focused 
on  determining  the  intensity  and  scale  lengths  required  to  accelerate  electrons  to  energies 
of  5  MeV  or  more.  It  was  found  in  earlier  work,  that  a  very  high  intensity  of  20  W/cm^ 
is  required  for  this  purpose  in  the  nighttime  ionosphere,  while  the  optimum  frequency  is 
roughly  twice  the  electron  gyrofrequency.  During  the  present  contract  period,  a  diffusion 
code  which  solves  the  Fokker-Planck  equation  along  the  Hamiltonian  surfaces  to  determine 
the  evolution  of  the  distribution  function  was  completed.  Intensities  roughly  a  factor  of 
four  over  threshold  are  required  to  yield  good  agreement  between  the  particle  and  diffusion 
codes.  The  diffusion  code  runs  far  more  rapidly  than  the  particle  code  and  allows  us  to 
make  broad  studies  of  parameter  space.  The  work  led  to  two  manuscripts  submitted  for 
publication  to  Physical  Review  Letters  and  Journal  of  Geophysical  Research.  These  are 
included  in  Appendix  A. 

The  work  on  the  Alfven  maser  started  by  determining  the  mode  structure  of  low- 
order  shear  Alfven  modes  in  the  magnetosphere,  including  their  transverse  structure,  and 
ensuring  that  their  coupling  to  the  compressive  modes  is  weak.  With  the  study  of  the 
mode  structure  completed,  resonant  acceleration  mechanisms  to  determine  if  the  Alfven 
wave  alone  is  sufficient  to  precipitate  high  energy  ions  were  investigated.  The  results 
indicate  that  localized,  low-orde:  MHD  shear  Alfven  waves  are  not  sufficient  to  precipitate 
a  large  number  of  particles.  It  was,  however,  found  that  parallel  electric  fields  caused  by 
plasma  kinetic  effects  in  conjunction  with  values  of  0.1  can  induce  significant  high 

energy  particle  precipitation  for  electric  fields  above  a  threshold  value  which  was  derived. 
Such  waves  can  propagate  inside  natural  or  induced  magnetospheric  ducts.  We  note  that 
the  physical  mechanism  of  resonant  excitation  described  here  is  quite  different  physically 
from  previous  schemes  discussed  by  Trakhtengerts  and  co-workers. 
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II.  IONOSPHERIC  DIFFUSION 


A.  Test-Particle  Results 

The  diffusion  code  which  we  have  written  is  based  on  the  concepts  of  stochristicity  and 
resonance  overlap.^  These  concepts  have  played  an  increasingly  important  role  in  plasma 
physics  during  the  last  decade.  Our  basic  particle  Hamiltonian  may  be  written^ 

H  =  [(cP  —  qA)^  -|-  +  5$,  (2-1) 


where 


k»  ,  ,  kj_ 

A  =  Ai  —  sintpex  +  A2  cos V’Cy  —  -4i  —  sin V’Sz  +  xBoCy, 

K  K 


$  =  $0  sin?/?,  (2-2) 

xjj  —  k^x  +  A:||2  —  tut, 

and  q  is  the  particle’s  charge.  This  Hamiltonian  represents  single  particle  in  an  applied 
electromagnetic  wave  propagating  obliquely  with  respect  to  a  magnetic  field.  We  now  use 
a  series  of  canoniced  transformations  to  reduce  the  Hamiltonian  to  the  form^ 


The  zero-order  Hamiltonian  is  given  by 


u _ 2^■,  ,  pi  ,  pI  w 

Ho  —  me  (  1  -I-  — — ■  H - j-y  1  -  — 


2  tu 

=  me  70  -  — p||, 
fcji 

where  p±_  and  p||  are  momentume  perpendicular  and  parallel  to  the  magnetic  field,  Eind  the 
first-order  Hamiltonian  is  given  by 


/f.=  ^  Zi  sin(A:||Z  -|-  W), 


where  0  is  the  gyroangle  and 


mc2  [  q  P\\  .  q  m ,  n 

- <  Cl  — — sma-— ■; — cosa  Ji{k_ip) 

7o  I  |?|  ^'xc 

+  ^2  —  Jl(k±p}  +  (3-^^'roJl{kxp)\, 
me  |g|  J 


0 


with 


Here,  a  is  the  angle  of  wave  propagation  with  respect  to  the  magnetic  field  and  Q  Is  the 
magnitude  of  the  nonrelativisitic  gyrofrequency. 

Resonances  occur  when  k±  -{■  Id  =  0  or 

(2,8) 

mjo  To 

a  result  which  is  well-known  from  linear  plasma  theory.  The  resonance  surfaces  are  elliptic 
when  the  iif-surfaces  axe  hyperbolic  and  vice-versa.  If  the  wave  amplitude  is  small,  than 
an  electron  will  be  affected  by  at  most  one  resonance.  In  this  case,  all  but  that  single 
resonance’s  contribution  to  the  Hamiltonian  can  be  ignored,  and  the  Hamiltonian  becomes 

H  =  Ho{pi\,p±)  -f-  Zi  sin(^|i2),  (2.9) 

where  r  has  been  re-defined  to  absorb  W.  In  this  re-defined  coordinate  system,  the  reso¬ 
nance  occurs  at  a  point  (piid  P±r)  where  i  =  0.  Writing 


^^o(P||,PJ.)  ^  ffoCPjr.Pir)  +  X 


2  dpi  ' 


(2.10) 


and  letting 


M-'  =  ^ 


we  find  that  the  Hamiltonian  becomes 


(2.11) 


:(Pil  -P||r)^  +  Zisin{k\iz), 


(2.12) 


whicii  is  just  the  standard  pendulum  Hamiltonian.*  The  trapping  width  and  the  bounce 
frequency  are  given  by 

Api,  =  4|MZ,|*/'  (2.13) 


016  =  A:|||Z,/M!*/=, 


(2.14) 


respectively. 


It  is  important  to  note  that  the  term  “trapped”  has  a  very  specific  meaning:  it  refers 
to  orbits  close  enough  to  resonance  .so  that  within  the  rest  frame  of  the  resonance  (not 
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the  wave)  the  electrons  undergo  libratory  rather  than  rotatory  mocion.  Electrons  iiot 
clos:  to  a  resonance  still  undergo  acceleration  but  it  is  far  weaker  than  the  acceleration 
undergone  by  those  which  are  trapped.  Standard  rf  heating  schemes  make  use  of  the 
resonant  acceleration  just  described.  Of  course,  acceleration  is  limited  by  the  trapping 
width.  While  this  limitation  is  of  no  importance  in  schemes  to  heat  a  bulk  plasma  to 
moderate  energies,  it  is  of  great  importance  in  schemes  to  heat  electrons  to  ultrarelativistic 
energies.  In  this  latter  case,  a  stochastic  mechanism  is  needed  to  break  the  constraint 
imposed  by  the  finite  trapping  widths  of  the  individual  resonances.  The  separation  on  an 
//-surface  between  two  adjacent  resonances  is  given  by 


<5p||  = 


(2.15) 


The  stochasticity  threshold  condition  based  on  the  overlap  criterion  now  becomes* 

> 


n 


1. 


(2.16) 


Once  this  threshold  has  been  exceeded,  it  is  possible  to  accelerate  electrons  from  resoncince 
to  resonance.  We  can  determine  in  this  way  the  intensity  required  to  accelerate  ionospheric 
electrons  to  10  MeV  or  more.  These  calculations  may  be  found  in  reference  2. 


B.  Diffusion  Equation 

One  can  obtain  the  basic  diffusion  equation  using  either  a  Fokker-Planck  approach^"® 
or  a  quasilinear  approach.®”^  In  either  case,  one  obtains  the  equation 


dt  dv^  dv^ 


(2.17) 


The  Hamiltonian  is  a  function  of  p|[  and  px;  thus,  it  is  convenient  to  steurt  with  the 
relativistic  diffusion  equation  in  cylindrical  coordinates 


df 


1  d 


dt  px  dpi.  r 
where 


c)  d 

(^1111  tit:  + 


(2.18) 


^!!ll  =  / 

Jo 

-  r 

Jo 

D\\l.  =  Dn 


dr  (piKOpiK^  +  t)), 

dr  {p i_{t)p i_{t  -h  t)}, 

=  /  dr(pj_(t)pii(t  +  r)), 
Jo 


(2.19) 
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The  Jini^ular  brackets  indicate  an  average  over  the  gyroangle.  It  is  assumed  that  iiny 
initially  inhvjinc>gpneity  in  the  angular  distribution  is  rapidly  smoothed  out.  Integratmg 
over  unperturbed  orbits  -  the  standard  approxiuiation  of  fi”'''-ilinear  theory — we  find 


1=  00 

/=  — cso 

=  i?x||  =  ^hmQ  Y  —ZfsL  -  ^  ■ 

2  P  I  V  7717 


(2.20) 


tiipii  in 


where  Zi  are  resonance  coefficients  from  Eq.  (2.6)  and  ^(x)  indicates  the  usual  Dirac  6- 
function. 

These  diffusion  coefficients  are  not  useful  as  they  stand  because  they  do  not  take  into 
account  the  decorrelation  due  to  s.iochasticity  which  results  in  their  broadening.®  Here,  we 
will  assume  that  the  ^-functions  are  broadened  by  an  amount  which  equals  the  trapping 
widths  of  the  undisturbed  resonances.  This  assumption  works  well  in  practice  as  will 
shortly  be  apparent  from  our  results.  Thus,  we  replace  the  ^-function  by  a  function  f(x), 
where 


if|x|<2u;4 
~\0.  if  ixi  >  2w4 


The  quantity  u-'j  corresponds  to  the  bounce  frequency  in  the  /the  resonance. 

Because  the  Hamiltonian  if  is  a  conserved  quantity,  the  diffusion  proceeds  along  a 
constant  ii-surface,  and  we  can  reduce  the  diffusion  equation  to  one  dimension.  We  replace 
the  coordinates  (pipPx)  with  (if  ,7)  and  obtain  the  diffusion  equation  in  the  form 


^  -  i  —  f  n 

dt  757  V 


where 


D  -- _ 


(2.22) 


(2.23) 


Xotp  that  while  the  problem  has  been  reduced  to  one  dimension,  it  is  parameterized  by 
H.  and  Eff.  (2.23)  must  be  solved  for  the  entire  range  of  ii-values  imjic^rtant  in  any  given 
I'lrfdjlern. 


W^^  now  must  discreu.ize  Eqs.  (2.22)  and  (2.23)  so  that  they  can  be  solved  numerically. 
The  v.-ay  in  wliich  we  have  defined  implies  that  its  derivative  is  everywhere  zero  except 


on  a  countable  set  of  points  where  it  becomes  infinite:  hence,  the  discretization  requires 
great  care.  We  use  a  finite  element  discretization  approach  which  effectively  averages  the 
derivatives,  including  the  A-functions,  over  a  range  in  7.^’'°  Multiplying  Eq.  (2.22)  by 
where  g{y)  is  a  weight  'unction  to  be  made  explicit  shortly,  and  integrating  over 
the  domain  T  =  (7min,7max),  we  find 

9 

=  (2-24) 

where  wi?  have  used  integration  by  parts  and  pcirticle  conservation  to  eliminate  the  .second 
derivative  of  /  with  respect  to  7.  Dividing  the  range  F  into  N  ecjual  sized  elements,  we 
may  write  /  as 

N 

f  =  (2.25) 

J=1 

where  the  K’jiy)  are  the  basis  functions  for  the  elements,  defined  so  that 

{(7-7t-i/(7i-T.-i,  7i-i  <  7  <  7. 

(7.+1  -  7)/(7.+i  -  7i).  7.  <7<7.+  i  (2.26) 

0,  otherwise 

Using  the  1/1,  as  weighting  functions  and  defining  the  matrices  A  and  B,  such  that 

=  =  (2.27) 

we  find 

A#* 

A  •  —  =  B  •  f,  (2.28) 

where  f  indicates  the  vector  consisting  of  the  fj.  Since  A  and  B  are  both  tridiagonal, 

Ecp  (2.27)  is  not  difficult  to  solve. 

We  liave  tested  the  code  by  running  it  for  cases  where  the  diffusion  coefficient  has 
analytical  solutions  and  found  the  agreement  between  numerical  and  analytical  results  to 
be  excellent.  We  have  also  found  the  number  of  particles  to  be  typically  constant  to  better 
than  one  part  in  10“^. 

C.  Computational  Results 

In  all  our  numerical  calculations,  we  have  taken  the  initial  distribution  to  be  a  (5-  » 

function  at  7  =  1  and  we  have  varied  our  parameters  over  a  wifle  range  of  wave  ener¬ 
gies.  The  parameter  values  we  have  chosen  are  consistent  with  the  daytime  and  nighttime 
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>  iM  iiii  1 1'^  ‘ '  I:;  i  r';r:a;  i;;^t  anct's.  we  liavf  coaipaicd  sni‘3;lc  pJLrtic'lc  raic'jlati(jnh  tlic 

a.-Jia  tac  dithiMtJU  I'o.ic.  In  these  eaicuiatioiis.  we  liave  used  25G  pmticles.  \V liile 
•iu-'i  ijt.i  »>!'  parta  ies  ;.s  not  laiare  e:ioui;;li  to  .u,!''--  us  a  good  solution  fen'  tlie  oe'erall  evo- 
intani  in  tia'  di.-'tr;init.ju  tunetions,  it  is  sufficient  to  allow  ns  to  solve  for  crucal  statistical 
nauanieters  such  .  '  the  average  distance  travelled  by  the  electrons  and  the  first  two  mo- 
:ne:its:  -•  ~  ttnd  '■  —  'vhere  70-  initial  7-value,  equals  one.  We  carry  out  our 

studu's  for  a  length  itf  time  \lt  —  30,000.  ■which  corresponds  typically  to  travel  distance  of 
alio'it  20  km  for  the  ele'c/trons." 

Firs',  '.ve  considi'r  tlie  diffusion  due  to  a  R-x  modi'  with  =  1.97.  o  =  S0°.  and 

il  —  0,3.  The  ^■alue  Q  =  0.3  correspemds  to  the  nighttime  ionosphere  at  130  km." 
Ti.e  fiire.siioid.  for  the  onsrr  of  stochasiicity  i.s  fthrs  0.2. 

W'e  '.'airy  out  (.mr  comparison  for  ave  amplitudes  ranging  from  0.19  to  10.0.  The 
•.■;ir;a':oii  (,if  several  quantities  is  shown  in  Fig.  1  for  f  =  0.19.  The  quantity  {(7  —  70)^) 

>  ^iTii  to  increase  n  an  approximately  linear  fashion  as  expected  from  the  Fokker-Planck 

'■'fp,  For  this  set  of  parameters.  eventually  becomes  very  small,  and  ((7  —  7o)^) 

tiattens  out.  but  that  does  not  occur  until  Qt  ^  3000.  The  quantity  D/Dqi,  which  is 
the  ratio  of  the  orbit  calculation  of  ((7  —  7o)^)  to  that  obtained  from  the  diffusion  code, 
deviates  significantly  from  unity  at  earlier  times  but  settles  down  close  to  its  asymptotic 
value  for  Ut  ^  750.  Tlie  d.eviation  from  quuoilinear  behc\'ior  at  early  times  is  largely  due 
to  two  effects,  stickiness  and  excess  coherence  due  to  a  small  number  of  resonances.^  The 
quantities  \z')  and  corrc.spond  to  the  average  distance  travelled  for  all  the  electrons 

and  only  those  electrons  with  7  =  10.  Tla'se  two  quantities  become  increasingly  equal  at 
large  times  as  an  increasing  fraction  of  the  electrons  surpass  the  \'alue  7  —  10.  We  can 
t.'.stimate  ttie  value  /f  to  ijc' 

m  =  (7  -  ~  ,  (2.29) 

wiiiaii  liold.^  for  highl}'  I'  lai.ivi.-ric  electrons.  The  solid  line  in  Figs.  Id-e  is  given  by 
Fq,  2.29’.  and  riie  de-.-iarion  indicates  that  not  all  the  electrons  have  liecome  relativistic 

L:.e  tiaeryiii  (it  iject :  ons  iia\’e  7'  K)  is  [ilotr.i'd  in  Fig.  2.  The  cnrvt's  obtained  from 
r.ne  di.'iuuon  eode  i  -(/iiii  line/  and  rhc'  jjarfa'le  code  (flotte<l  line)  agree  cjuite  well.  Similar 
leadns  are  -iiown  m  rigs.  3  0  for  wave  tinqilitudes  given  by  f  =  3.99  and  f  ~  10.0.  The 
agreement  only  iinnroves  in  these  <...uses.  lii  Fig.  7,  we  show  Dj Dqi  in  the  range  f  =  0.19 
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to  e  =  10.0.  The  quasilineai-  tlieory  embodied  in  our  diffusion  code  evidently  works  well 
over  rhi.-;  rane,c. 

W’e  next  consider  the  diffusion  due  to  the  L-o  mode  with  —  2.G,  a  =  S0°,  and 

~pi\l  0.3.  The  threshold  for  stochastic  acceleration  is  given  roughly  by  Cthrs  =  0.4. 
Results  for  several  different  energies  can  be  found  in  Figs.  8-10.  In  this  case,  a  substantial 
fraction  of  the  electrons  are  not  stochastic  even  when  the  threshold  is  exceeded,  and  that 
has  an  important  effect  on  Dj Dqi,  (z),  and  {z)-y>io-  Nonetheless,  the  fraction  of  electrons 
with  y  :>  10  agrees  well  with  both  methods. 

The  results  of  our  study  of  electron  acceleration  generrlly  confirm  the  results  obtained 
unng  simpler  methods  in  reference  2.  In  particular,  large  intensities  of  roughly  20  W/cm* 
are  needed  to  reach  the  required  stochasticity  thresholds  and  acceleration  distance  on  the 
order  of  20  km  are  leeded  in  order  to  accelerate  a  substantial  fraction  of  electrons  beyond 
5  MeV.  1  he  power  requirements  appear  to  be  well  beyond  what  can  be  managed  with 
present-day  devices.  Other  mechanisms  based  on  field  inhomogeneities,  pre-acceleration, 
and  precipitation  of  already  high  energy  particles  remain  to  be  more  thoroughly  exp'ored. 
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III.  ALFVEN  MASER 


A.  Overview 

Tlie  Alfven  r.inser  i  A-iuaser),  as  first  described  by  Trakhtengerts  and  co-vvorkers^^”^^. 
is  based  on  enlirntcing  the  reflectivity  of  the  Alfven  vvav^e  within  a  narrow  tube  (  ~  10  km 
X  10  km)  along  the  field  lines  in  the  magnetosphere.  Tliis  change  in  reflectivity  could 
be  accomplished  by  rf  heating  of  the  ionosphere  in  the  region  near  the  poles  where  the 
magnetic  field  lines  are  tied  down.  Trakhtengerts  and  his  colleagues  focused  primarily  on 
metliocis  for  generating  the  ion  cyclotron  instability,  as  was  approjiriate  for  the  problem 
of  precipitating  relati\ely  low  energy  ions  which  he  was  considering.  Our  focus  has  been 
fhe  precipitation  of  higher  energy  ions  and  the  exploration  of  resonant  methods  which  can 
accomplish  thi-s  prccijiitation.  Papadopoulos.  ct  have  shown  that  the  mechanisms 

discu.>se(i  by  Trakhtengerts  and  co-workers  are  insufficient  for  this  purpose  for  several 
reasons,  inunely: 

fi)  Xo  specific  wave  mode  that  interacts  strongly  with  energetic  protons  (f  >  1 
MeV)  was  identified. 

(ii  )  The  models  rely  on  (juasilinear  wavepacket  theory.  As  a  result  their  applica¬ 
bility  is  limited  to  proton  energies  below  1-2  MeV. 

iiii)  The  transverse  variation  of  the  medium  and  finite  ion  gyroradius  effects  were 
ignored. 

.As  a  result,  we  concentrated  on  scenarios  which  excite  the  highly  localized  (rn  — +  oo) 
poloidal  oscillations  of  the  magnetospheric  cavity  at  the  appropriate  L  value  by  modulated 
RF  heating  of  the  F-region  ionosphere.  These  modes  are  thought  to  be  excited  naturally 
by  inverted  ion  distribution  functions,*^  *®  ehtetron  beams,'"  surface  waves  due  to  the 
Kelvin-Helmholtz  instal'ility,'®''^  and/or  external  impulses.""'^'  These  are  monochromatic 
fiehl  line  oscillations  and  have  often  been  observed  in  conjunction  with  energetic  particle 
nreci[Mtation.  We  have  named  the  process  Field  Resonance'  Induced  Precipitation  (FRIP). 
.X'jtice  rliat  riie  pi.y.-ax  o.'  f'PIP  is  completely  different  than  the  physics  of  the  Trakht- 
entmits  .\-.\Ia..ser.  lia.-rd  nn  onasilinear  diffusion,  both,  in  the  modi'  excitation  j^roci'sses  as 
well  as  the  part  ich’  pa  eci  p: '  a  t  mn  pi  ocess.  The  difleri'iices  will  be  einiihasizefl  hater  on.  The 
two  important  issues  related  to  the  feasibility  of  FRIP  are 
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(i)  Identification  and  verification  of  the  physics  related  to  the  induced  energetic 
proton  precipitation. 

(ii)  The  requirements  and  location  for  exciting  the  appropriate  field  line  reso¬ 
nances 

These  are  discussed  in  the  next  section. 

B.  Radiation  Belt  Dynamics 

Before  entering  the  discussion  of  the  physics  issued  relat^^d  to  the  FRIP  concept  we 
briefly  present  the  parameters  that  control  the  trapping  and  the  precipitation  of  the  ener¬ 
getic  protons  in  the  earth’s  radiation  belts.  Figure  11  shows  the  trajectory  of  a  trapped 
particle  in  the  eaxth’s  dipole  field.  Particles  trapped  in  the  earth’s  RB  execute  three  types 
of  motions.  First,  a  gyration  about  <’he  magnetic  field  with  the  local  gyrofrequency  Clc, 
second  a  longitudinal  motion  governed  by  the  effective  potential  energy, 


U||  =  nB{s) 

where  B  is  the  magnetic  field,  s  is  the  coordinate  along  a  line  of  force  and 


(3.1) 


2  2-2 
V  sin  a 

W  ^  2B 


(3.2) 


is  the  orbital  magnetic  moment;  a  is  the  angle  between  the  velocity  v  and  the  field  B.  This 
motion  results  in  the  particle  bouncing  between  mirror  points  along  the  field  line  with  a 
frequency  uJt,.  Third  is  a  slow  azimuthal  drift.  These  motions  are  shown  in  Fig.  11.  The 
value  of  the  magnetic  field  Be  at  the  equator  is  given  approximately  by 


BeiL)  =  Y3  ’ 

while  the  average  value  ot  the  gyrofrequency  in  the  dipole  field  is  given  by 


(3.3) 


(I,  -  n 


cE 


F(a) 


(3.4) 
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where  Q^E  =  3  x  sec~^ ,  G(a)  =  1.3  — .56  sin  a,  and  F(a)  ~  (—0.2  sin^  a  +  0.7 sin  a  + 

2)/  sin'  a.  In  the  absence  of  interactions,  the  magnetic  moment  //  is  an  adiabatic  invariant. 
The  same  holds  true  for  the  action 


J  —  j)  V  ■  ds  , 

which  is  associated  with  the  bounce  motion.  There  is  also  a  third  adiabatic  invariant,  the 
flux  invariant  associated  vvith  the  drift  motion. 

C.  Induced  Precipitation 

Energetic  protons  will  be  confined  in  the  radiation  belts  as  long  as  their  adiabatic 
invariants  are  conserved.  The  rate  at  which  protons  enter  the  loss  cone  and  precipitate 
depends  mainly  on  the  time  scale  on  which  the  adiabatic  invariant  J  cheinges.  Two  general 
types  of  processes  exist.  First,  low  frequency,  random  fluctuations  will  lead  to  a  quasilinear, 
diffusive  random  walk  in  J .  Second,  coherent  fields  can  lead  to  resonant  excitation  and 
stochasticity.  It  is  this  latter  type  of  process  which  we  are  proposing  to  use,  rather  than 
the  former  type  proposed  by  Trakhtergents  and  co-workers.^^’^^ 

The  guiding  center  of  a  trapped  particle  oscillates  between  its  mirror  points  with  a 
frequency  given  approximately  by 


u,  =  2.2 


y/T(MeV) 


(3.5) 


In  general,  the  magnetic  field  does  not  vary  quadratically  v/ith  distance  along  the  field 
line,  leading  to  a  significant  anharmonicity  in  the  bounce  frequency.  In  other  words, 
where  3£  is  the  particle’s  pitch  angle  at  the  equator,  for  particles  with  the 
same  value  of  T.  However,  this  anharmonicity  does  not  significantly  alter  the  results  and 
does  complicate  the  discussion.  So,  we  ignore  it  here.  Our  equation  of  motion  is  given  by 


S  +  CJj  5  =  0  , 


which  comes  from  the  Hamiltonian 


1 


H„  =  ^  =  0 

Irn  1 


(3.6) 


(3.7) 
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Making  the  action-angle  transformation. 


cos^  ,  s  =  sin6)  . 


(3.S) 


the  Hamiltonian  becomes 


Ho  =  ujbl 


(3.9) 


with  equations  of  motion 


7  =  0,  6  =uji 


(3.10) 


The  presence  of  electromagnetic  forces  with  components  pareJlel  to  B.  will  change 
the  value  of  J .  Forces  parallel  to  B  can  be  due  to  parallel  electric  fields  E\\  or  due  to 
compressional  magnetic  perturbations  6|| .  The  latter  represent  a  form  of  Fermi  acceleration. 
In  the  presence  of  such  forces,  the  equation  of  motion  becomes 


..  ,  2  ef;||(s,t)  56||(s,0 

s  ^ - p - - - 

m  os 


(3.11) 


The  strongest  interaction  will  of  course  occur  for  forces  oscillating  near 

We  now  suppose  that  the  fields  are  fluctuating  with  random  time  variations.  Such 
fluctuations  induce  random  variations  in  the  oscillation  amplitude,  causing  a  random  walk 
of  a  particle’s  J  value  and  a  redistribution  of  J  described  by  a  Fokker-Planck  equation. 
Ro’oertc  and  Schultz^^  showed  that  spatial  variations  allow  the  oscillator  to  interact  with 
the  fields  at  harmonics  of  the  bounce  frequency.  However,  the  overall  effectiveness  of  the 
bounce  resonance  processes  is  reduced  for  oscillation  amplitudes  large  with  respect  to  the 
scale  of  the  spatial  variations.  For  situations  where  p  =  const,  the  Fokker-Planck  equation 
describing  the  evolution  of  the  distribution  function  /(e|i,ex)  is 


^/(g)  ^  ^  p  ^  r 

dt  den  den  ^ 


(3.12a) 


_  2  ^  2  [°°  t2/i.  /Vp(A:||,nU.'),) 

,>  =  1  ^11 


(3.126) 
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where  is  the  maximum  value  of  5  in  the  unperturbed  oscillation,  the  are  Bessel 
functions,  and  Kp  is  the  relevant  power  spectrum.  These  relations  can  be  derived  from 
Ecj.  (3.11)  by  using  second  order  perturbation  theory  and  defining  D  as 


1  <  (Ae||)^  > 

2  T 


(3.13) 


Of  greater  interest  to  us  is  the  case  where  we  may  assume  £'||(s,t)  =  Eg  cos{ks  —  Lut). 
In  later  sections,  we  will  explicitly  consider  Alfven  waves  in  which  case 


ks I  — ’  (3-14) 

Jo 

and  Eo  becomes  a  weakly  varying  function  of  s.  The  effect  is  to  add  harmonics  of  the 
fundamental  wavenumber 


Eo  cos(^'s  —  wf)  — >  ^  cos  ((2n  +  l)^:s  —  u;f]  ,  (3.15) 

n=0 

which  has  little  impact  on  the  final  result.  Our  equation  of  motion  now  becomes 


2  ,,  . 

- cos{ks  —  ujt)  , 

m 


(3.16) 


Letting  a  =  ks  and  r  =  Eq.  (3.16)  becomes 


a  +  a  = 


eEgk 


mu)i 


cos(cr  —  Or)  , 


(3.17) 


where  0  =  u;  jub-  This  equation  of  motion  can  be  derived  from  the  Hamiltonian 


H  +  ^  +  9.R  -  sin(cr  -  <^) 


(3.18) 


.vhere  (TZ.cp)  are  a  canonical  pair  of  variables.  Making  the  action-angle  transformation 


^  =  (2/)’/%in^  ,  =  (27)'/^cos^  , 


(3.19) 
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so  that  in  physical  units  I  =  (U  —  where  U  is  the  kinetic  energy  and  Dg  is 

the  magnetic  field  at  the  equation,  our  Hamiltonian  now  becomes 


H  =  I  +  Q.R  —  e  sin 


(2/)‘/'''sin^  -  4> 


(3.20) 


where  e  =  eEok/mu;'^  which  is,  to  within  a  phase,  the  well-known  Hamiltonian  considered 
by  Karney"^  in  the  context  of  rf  heating.  For  a  fixed  value  of  e  and  Q,  Karney  showed  that 
the  orbits  are  stochastic  in  the  range 


,  (3.21) 

where  the  latter  value  holds  exactly  in  the  limit  (2iy^^  »  Q,  but  yields  a  reasonable 
approximation  even  when  {2iy^^  ~  Q.  An  illustration  of  the  stochastic  region  for  a  fixed 
value  of  fl  is  shown  in  Fig.  12. 

As  an  example,  we  find  that  if  Q  =  10,  then  the  threshold  for  stochasticity  is  approx¬ 
imately  at  e  =  1;  however,  larger  values  are  required  to  obteiin  significant  acceleration.  In 
physical  units,  our  conditions  become 


U! 

iJB 


TnuJBVAE  <-^B  J  \^/  V  ^^BVaE  ^B  / 


(3.22) 


where  Vae  =  va  at  the  equator.  Writing 


eEp 

mUBVAE 

2T 

3mv\E 


-  n  nwgo  r5/2 

yr(Mev) 

-  3.3  [T(Mev)]  , 


(3.23) 


we  see  that  at  L  =  4,  and  setting  as  before  12  =  5,  a  field  of  7  inV/m  is  needed  to  achieve 
f  =  5  for  a  200  kev  proton.  .A.t  this  value  of  E,  protons  are  stochastic  if  30  <  /  <  430 
and  our  200  kev  proton  has  /  =  42.  We  see  that  if  e  =  5  and  12  =  10  are  maintained  by 
ramping  up  the  E-field  and  the  frequency,  a  parallel  energy  gain  of  10  is  possible,  leading 
to  precipitation  of  most  200  kev  protons  One  can  think  of  many  other  possible  scenarios 
for  precipitation.  .A.s  will  be  discussed  later,  such  fields  can  be  excited  due  to  finite  ion 
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gyroraclius  effects  when  Alfven  field  aligned  resonances  are  driven  by  periodic  heating  of 
the  F-region. 


D.  Geomagnetic  Field  Tube  Oscillations 

A  number  of  authors"^”^®  have  calculated  the  eigenfrequencies  of  magnetic  field  lines 
for  various  distributions  of  u  4  along  the  magnetic  field  lines  using  the  boundary  condition 
Eis  =  ±1)  =  0.  This  boundary  condition  is  correct  only  for  sufficiently  large  ionospheric 
conductivities.  The  resonant  frequencies  are  independent  of  the  conductivity  only  when 
v,\(s)  =  const  and  the  conjugate  ionospheres  are  symmetric.  In  this  case,  the  fundamental 
resonance  frequency  uJoiL)  is  given  by 


XV  A 


(3.24) 


where  Tq  is  the  earth  radius.  A  complete  analysis  using  the  inhomogeneous  wave  equation 
was  performed  by  Arykov  and  Maltsev. The  wave  equation  for  guided  Alfven  waves  is 
given  by 


^ _ ^  _  1  /  E  \ 

ds'^  s/Wo  v\{s)dt'^  \y/Sl) 


(3.25) 


In  Eq.  (3.25),  the  field  curvature  and  the  transverse  dependence  of  va  are  neglected.  This 
is  correct  for  relatively  small  transverse  dimensions.^*  The  resonator  properties  can  be 
described  by  introducing  the  function  A,  defined  as^^ 


AE  ^  [b;  AE 

Eo  V  Bo.  2Epn  ’ 


(3.26) 


where  AE  is  the  change  of  the  ambient  electric  field  Eg  in  the  ionosphere  if  the  Pedersen 
conductivity  varies  from  by  AE  in  the  northern  hemisphere;  Boi  is  the  magnetic 

field  on  tiie  ionosphere.  Note  that  if  A  =  1,  then  AE  equals  the  polarization  field  in  the 
modified  area,  on  condition  that  there  is  vacuum  above  the  ionosphere.  .Assuming  that 
AE  varies  periodically  with  frequency  lv,  the  equation  describing  the  excitation  factor  .4 
A 


d'^  up- 

- .4  4 - 

c>s2  ^  vA^ 


.4  =  0 


(3.27) 
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with  the  boundary  conditions 


.-1  + 


dA 


47ria,’Sps  ds 


0  s  =  -1  , 


dA  , 

A  +  - ;7—  =0  s  —  +l  . 


(3.2Sa) 


(3.286) 


47rtcjSp/v  ds 

The  function  ^4  depends  only  on  the  frequency  and  on  the  characteristics  of  the  resonator 
only.  Equation  (3.28a)  gives  passive  reflection  of  the  guided  waves  from  the  southern 
conjugate,  while  (3.28b)  describes  reflection  and  generation.  For  the  magnetosphere,  an 
analogous  excitation  factor  Ab  defined  by 


AB  =  -A, 


c  1  Bo  X  SqAE 


Va  yJBoBoi 

can  be  introduced.  Following  Arykov  and  Maltsev, we  approximate  va{s)  by 

=  vae  (^1  + 

In  this  case  Eq.  3.27  can  be  solved  to  give 


(3.29) 


A  = 


P  +  T]_  1  — 


where 


^3.30) 


(3.31) 


N,S 

n  -  _2±_ 

Rn,S  —  ^5 

V- 


^niuiEpisf^s 


I  ±  ia^/l  +  a^uR/vYE 


(3.32a) 


(3.326) 


1  +  a'^u;2/t;  j  t;  arctan(s/a)  , 


(3.32c) 


Uo  —  u(s  =  1) 


(3.32d) 
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In  order  to  analyze  the  situation  further  we  note  that  using  WKB  theory,  we  find 


TT  t^Vae 


2a  arctan(//a) 

J-l  tM(j) 


(3.33) 


where  u.'o  is  the  frequency  of  the  fundamental  mode.  We  now  introduce  the  inhomogeneity 
parameter  X  =  l/a  and  the  ratio  of  the  ionospheric  conductivity  to  the  wave  conductivity  at 
the  equator,  a  =  Dp/Su^,,  where  11,^;  =  jAirv^E-  If  A  =  0,  corresponding  to  a  homogeneous 
magnetosphere,  we  find 


A  = 


exp 

1 

. 

f 

+  Rs 

exp 

*(7  +  3)' 

1  -  Rt^Rsexp  1 

f27rt^'] 

1 

(3.34) 


w 


here 


-R/V,s  =  (1  -  -r  <’’iV,s)  •  (3.35) 

Aryhcv  and  Maltsev^^  calculate  the  value  of  A  in  the  conjugate  ionospheres.  Our  interest 
here  lies  on  the  equatorial  pleme  s  =  0.  The  results  are  shown  in  Fig.  13,  for  symmetric 
and  antisymmetric  ionospheres.  The  most  interesting  result  is  the  one  shown  in  Fig.  14 
which  shows  the  amplitude-frequency  chaxaeteristic  for  Ab  =  i{vA/i^){dA/ds).  A  strong 
resonance  is  not  possible  on  the  odd  harmonics,  but  only  on  the  even  ones,  i.e.,  u!  =  2nu}o- 

Using  the  simple  model,  ==  t;^o(ro/r)^/^,  where  =4.4x10®  cm/sec  is  the  Alfven 
speed  near  the  Earth’s  surface  and  is  the  Earth’s  radius,  we  find  from  integration  of 
Eq.  (3.33)  that  tJo  =  6-8  x  10"^  sec“^  when  X  =  4.  Hence,  in  the  example  of  the  previous 
section  in  which  iu  =  IOu^a  =  2.5  sec~^  we  find  that  n  ~  18. 

E.  Stability  of  the  Shear  Alfven  Wave 

In  this  section,  we  will  show  that  the  transverse  mode  structure  can  be  determined 
separately  from  the  longitudinal  structure  and  that  the  shear  Alfven  wave  is  stable.  We 
have  investigated  both  diffraction  and  mode-mode  coupling  and  shown  that  they  are  com¬ 
pletely  negligible  at  the  low  parallel  mode  numbers  of  the  shear  Alfven  wave.  Indeed,  the 
only  limitation  on  the  transverse  size  of  the  mode  is  that  it  must  be  large  compared  to  an 
ion  gyroradius  for  the  basic  dispersion  relation  to  be  valid. 
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Our  starting  point  is  the  equations 


V  X  B  =  -  — K  -E, 

c 

—  ^  —to,’  „ 

V  X  E  - B, 

c 


(3.36) 


where,  assuming  we  have  low  frequency  modes,  but  keeping  liigher-order  cyclotron  correc¬ 
tions,  we  have 


K  = 


0 


(3.37) 


where,  writing  E  =  (Ei,  E2,  ■E'3)^  •E'l  and  E2  are  transverse  to  the  zero-order  magnetic 
field  and  E3  is  directed  along  it.  At  low  frequencies,  E3  is  essentially  zero  for  consistency. 
To  drive  an  appropriate  Fresnel  criterion  for  Alfven  waves  in  magnetospheric  flux  tubes, 
we  first  ignore  geometric  effects  and  obtain 


To  lowest  order  in  /fif,  we  obtain  the  dispersion  relation 


(3.38) 


V 


2 

A 


/  >2  ,  ,2  ,  ,2 

U/  UJ  (jJ 

If 


when  '>  A:||  as  in  the  case  in  a  confined  mode.  We  then  find  that 


(3.39) 


d'.j  (x>  ^11 


(3.40) 


Thus,  the  more  we  reduce  the  transverse  dimensions  of  the  wave,  the  less  it  diffracts! 
We  conclude  that  diffraction  imposes  no  limit  on  the  transverse  width  of  the  mode.  The 
width  need  only  be  large  enough  to  validate  the  dispersion  relation,  i.e.,  it  must  be  large 
compared  to  the  ion  gyroradius.  This  radius  is  roughly  10  m  when  L  =  4.  Hence,  a 
transverse  width  of  a  kilometer  or  more  should  be  sufficient  to  provide  confinement. 
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(3.41) 


To  investigate  tlie  effect  of  mode  coupling,  we  first  note  that 

2 

=  ^nNkTD^  ~  10-* 

when  L  —  4.  Hence,  acoustic  coupling  is  completely  negligible.  To  consider  coupling  to 
the  magnetosonic  wave,  we  write  Eq.  (3.36)  in  dipole  (Radoski)  coordinates.*®  It  becomes 

=  —  (  — T~^‘' - )  “^ - (  ^  — T  )  ^‘t>'  (3.42a) 
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^  c 

1  +  (3.426) 


2 

(jJ  C 


02  „2 


1  ^  i_  rr  n 

T  7  huhju  = 

hi,hn  a iJ.  c 

1  d  iu 

dll  c 

where  for  any  quantity,  we  have  let 

A’(/i,  n,  <j))  =  X(/2)exp(tsi/  +■  im(t>). 

Defining  now, 


„  m  ^  s  ^ 

R\  =  7~Rv  — 

h,j,  h„ 

„  m  „  s  ^ 

h,p  h„ 


(3.42c) 

(3.42d) 


(3.43) 


(3.44) 


so  that  R\  corresponds  to  the  magnetosonic  wave  and  corresponds  to  the  shear  Alfven 
wave,  we  find  that 


J-A_LA 

dll  h^i  dll 


hu  fh4>dK\  1  d  -  s^hD 

^  hi,  diih,*,)  h\dii{m^hl  +  s'^hy) 
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■  i  H  ^  T  ‘ 

h,j,  d  h„\  2  d  msh„h,i, 


h^Ri 


hi,  dll  hij,  J  h'j.  dll  (s2/i2  -p  rn2/i2) 


h  ^i  R.‘2  , 


(3.45) 
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where  wc  neglect  the  vacuum  and  cyclotron  contributions  to  Eq.  (3.42).  If  either  rn  or  5 
is  equal  to  zero,  the  magnetosonic  wave  is  not  coupled  to  the  shear  Alfven  wave.  That  is 
because  curvature  induces  a  field  pointed  in  the  radial  direction  to  produce  a  component 
in  the  parallel  direction  while  a  field  pointed  in  the  azimuthal  direction  is  not  affected  by 
curvature.  Coupling  cannot  be  induced  in  either  case.  Noting  that  at  low  parallel  mode 
numbers  all  parallel  derivatives  d/h^d^  scale  as  1/ro,  and  assuming  as  well  that  A'xCq  ^  1. 
where 


2 


(3.46) 


we  find  that 


Rx  (3.47) 

Rd 

where  roughly  equals  the  transverse  width  of  the  confined  mode.  Thus,  the  intensity 
of  the  magnetosonic  wave  is  a  completely  negligible  fraction  of  the  shear  Alfven  wave  and 
leads  to  negligible  attenuation. 


——r 


(3.48) 


Thus,  we  may  set  Ri  =  0  and  all  dependence  on  m  and  s  disappears  in  Eq.  (3.42).  We 
may  in  principle  determine  the  entire  mode  structure  by  solving  one  pair  of  equations,  for 
instance 


1  d 

ht,h^  dfj. 


1  d 


and  determining  E^,  and  through  the  relations 


(3.49) 


(3.50) 
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These  results  imply  that  the  longitudinal  and  transverse  parts  separate  and  that  the 
transverse  portion  is  substantially  but  not  completely  arbitrary.  From  Eq.  (3.49),  we  may 
write  a  solution  which  satisfies  a  particular  set  of  boundary  conditions  as 


E4.  = 


(3.51) 


The  two  functions  gliji)  and  are  determined  from.  Eq.  13.49),  which  is  purely  lon¬ 

gitudinal.  Note  that  the  transverse  shape  factor  for  and  is  necessarily  the  same. 
Similarly,  we  may  write 


Et,  =  o), 

0). 


(3.52) 


Our  task  now  is  to  determine  what  constraints  Eq.  (3.50)  imposes  on  the  relationship 
between  fa  and  fj.  Using  Fourier  transform  methods,  we  can  how  that  this  equation 
implies 


(3.53) 

and.  what  is  more,  Eq.  (3.53)  implies  Eq.  (3.50).  Thus,  we  choose  either  fa  or  //j  arbitrarily, 
but  the  other  is  then  fixed  to  within  a  constant. 

.■\s  an  example,  we  may  consider 


fa(  V.  0)  =  /o  exp[-a(i/  -  i/o)‘  -  1)(0  -  0o)^] .  (3.54) 

whicii  corresponds  fi>  a  localizetl  .A.lfven  mode.  It  follows  that 

/j  =  0),  (3.55) 

a 

-o  that  the  relative  strength  of  these  *wo  shape  factors  is  determined  by  the  falloT  of  the 
overall  amplitude.  Equation  (3. .55)  can  be  viewer!  as  a  result  of  qua.si-neutrality.  In  this 
'■xarnrtie.  the  wave  is  linearly  jiolai'zed,  but  other  polarizations  result  when  we  allow  a  and 
h  to  be  roinplex.  Ill  particular,  if  we  set 
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a 


i  3.5G) 


—  i 


V2 


■A. 


b  =  i^.4. 
^/2 


\vc  find  that  the  resultinc;  wave  io  circularly  polarized.  As  a  final  j)oiat.  we  note  that  a 
circularly  syrnmetoric  shape  cannot  in  general  be  maintained  along  an  entire  Alfven  tube. 
To  obtain  a  circularly  symmetric  shape  at  a  single  point,  we  demand  that 


fc  ^  fj  -fl  [a(i'  -  I'oF  +  b{(p  -  (3.57) 

where  a  =  /i*^  and  b  =  at  that  point.  As  and  change  along  the  field  line,  the 
shape  will  become  assymmetric. 

We  inay  thus  determine  the  resonant  frequencies  by  solving  the  one- dimensional  lon¬ 
gitudinal  equations.  In  the  WKB  approximation,  we  obtain  the  frequencies  derived  in  the 
previous  section. 

F.  Kinetic  Alfven  Wave  and  Proton  Precipitation 

The  analysis  presented  in  the  previous  two  sections  was  based  on  MHD  equations, 
In  the  presence  of  transverse  gradients,  the  Alfven  mode  has  a  singularity.^*  *®  For  time 
varying  waves  this  effect  is  manifested  in  an  increase  in  the  transverse  wave  vector  of  the 
wave  described  by 


(9k  _ 

'm  "  '  It  ^ 

Under  magnetospheric  conditions  this  increase  in  the  wave-number  occurs  over  a  very 
fast  scale.’®  so  i,hat  the  transverse  wave-vector  /:x  becomes  larger  than  A'||.  In  order  to 
account  for  the  transverse  structure,  non-ideal  MHD  effects  must  l>e  incorporated.  nese 
are  the  electron  inertia  (c/a>)  and  the  ion  gyroradius  /?,.  Incorporating  these  effects,  the 
dispcvsiori  relation  becorne.s®* 


—  A:|l  1^4(1  .  (3.59) 

The  value  of  .1  depend.s  on  the  ion  temperature.  For  cold  ions.  J  <  and  we  find 

while  in  tlu-  opposite  limit  we  find  ,4  ~  In  the  latter  ca.se,  a  kinetic  analysis 

yields" 


99 


(3.60) 


j.2^^2 _ _ 

1  - /o(  A,)exp(-A,) 


+  Aj, 


where  A,  =  ,  Aj  =  \,T,ITe_,  and  is  the  modified  Bessel  function  of  the  first  kind. 

This  wave  is  known  as  the  kinetic  Alfven  wave  (KAW).  The  presence  of  a  weak  transverse 
dispersion  produces  a  small  group  velocity  across  B,  given  by 


(3.61) 

For  the  presence  of  weak  dispersion  can  trap  the  KAW  near  the  mini¬ 

mum  of  the  transverse  profile  of  Va  and  produce  an  Alfven  waveguide  for  the  KAW  similar 
to  the  one  discussed  in  the  previous  sections  for  the  shear  Alfven  wave.  Such  a  waveguide 
can  be  produced  by  the  plasmapause  as  well  as  by  magnetospheric  ducts.  Notice  that  Pci 
pulsations  have  been  interpreted  as  eigenmodes  of  such  waveguides. 

The  presence  of  guided  KAW  modes  allowed  us  to  identify  a  new  process  that  can 
break  the  second  adiabatic  invariance  of  multi-Mev  protons.  The  process  relies  on  the  faet 
that  the  KAW  has  an  electric  field  coxTiponent  Ez  parallel  to  the  magnetic  field.  It  is  given 

by3^ 


du> 

dTY 


—  iij  k 


u:SB± 

{Te/T,) 

k^  X  Cj 

1  + 

1  -  I„{klRf)e->‘l^i 

(Te/r.) 

(3.62) 


For  ik_iR,Y  1,  Te  <  T,,  we  find 


E\\  =  JklR^,)  .  (3.63) 

,4ssuming,  k^_R,  =  0.1  which  is  a  reasonable  value  for  the  Pci  pulsations,  we  find  that 
Ej_  ~  10  mV/cm  in  the  example  of  section  3B.  This  field  strength  can  be  obtained  from 
ground-based  transmitters. 
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IV.  EXPERIMENTAL  REQUIREMENTS 
A.  Concept  Summary 

The  TRIP  concept  (Fig.  15)  relies  on  injection  cf  large  amplitude  ULF  waves,  gen¬ 
erated  by  modulated  heating  of  the  ionosphere,  into  the  geomagnetic  force  tube  bounded 
by  the  modified  ionospheric  area.  If  the  modulation  frequency  equals  one  of  the  field  tube 
resonances  or  its  harmonic,  the  plasma  filled  force  tube  behaves  as  a  high  Q  cavity.  As  a 
result  the  excitation  level  of  its  resonant  modes  (i.e.,  sheEir  Alfven  waves)  increases  sub¬ 
stantially.  The  guided  Alfven  waves  are  incompressive  and  their  electric  field  is  transverse 
to  the  ambient  magnetic  field.  However,  magnetic  field  and  plasma  inhomogeneities,  finite 
ion  Larmor  radius  effects  and  dynamics  couple  the  shear  Alfven  wave,  to  a  kinetic  Alfven 
wave  (KAW).  The  kinetic  Alfven  wave  is  also  a  confined  mode  in  the  presence  of  plasma 
ducts. 

The  KAW  hats  an  electric  field  component  parallel  to  the  magnetic  field.  When  the 
value  of  the  parallel  electric  field  in  the  equatorial  region  exceeds  values  of  mV/m  protons 
with  hlev  energy  become  strochastic  (i.e.,  their  second  adiabatic  invariant  is  broken)  amd 
can  thus  precipitate. 

The  requirements  for  breakdown  of  the  second  adiabatic  invariant  and  induced  proton 
precipitation  can  be  approximated  by 


Wf,  =  2nuo 


(1) 


Eo  >  ~<jjILRe  (2) 

where  Eo  is  the  parallel  electric  field  of  the  KAW,  is  the  bounce  frequency  of  the 
precipitated  protons,  uJo  is  the  fundamental  frequency  of  Alfvenic  oscillations  of  the  field 
line  with  fixed  ends,  n  the  harmonic  number,  the  earth  radius  and  L  the  magnetic  shell 
number.  If  we  assume  a  simple  model  of  the  equatorial  values  in  the  radiation  belts 

B  =  .3j^G  (3) 

n.,  =  5  X  10'‘^#/cm^  (4) 
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and  take 


u,’6  « 


(5) 


We  can  express  (1)  and  (2)  in  a  quantitative  fashion.  The  value  of  n  required  for  precipi¬ 
tating  protons  of  energy  T(MeV)  at  a  value  L  is  given  by 


L  I  T 
'^"S.SVMeV 

and  the  associated  value  of  Eg  is  given  by 

Figure  16  gives  the  value  of  fo  =  tao/27r  as  a  function  of  L  for  the  model  given  by  (3)  and 

(4). 

For  example,  in  order  to  precipitate  1  Mev  protons  at  L  =  6  we  should  drive  the  shell 
at  its  fundamental  frequency  /<,  =  .03  Hz,  with  sufficient  strength  that  the  value  of  Eo 
at  the  equator  exceeds  2  mV/m.  To  precipitate  4  Mev  protons  we  should  drive  the  field 
at  its  second  harmonic,  i.e.,  .06  Hz,  and  the  required  value  of  Eg  would  be  8  mV/m,  etc. 
The  practical  implementation  of  the  system  depends  on  the  efficiency  with  whi-'h  ground 
based  HF  power  can  be  converted  into  KAW  power  in  the  magnetospheric  cavity.  This  is 
examined  next. 


B.  System  Requirements 

In  terms  of  implementing  FRIP,  the  major  issue  concerns  the  efficiency  with  which 
modulated  HF  power  injected  in  the  ionosphere  is  transformed  into  power  in  the  KAW 
in  Tie  equatorial  region.  We  present  below  estimates  based  on  modifying  the  ionospheric 
conductivity  of  the  lower  ionosphere  (D-E  region)  (Figure  15). 

For  a  dipole  magnetic  field  the  value  of  AE  of  the  electric  field  in  the  equatorial  plane 
i.s  given  by 


AE  -  A 


1  1 
cos^  ip 


AS 

V 
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where  Eg  is  the  ionospheric  field,  cp  the  latitude  and  A  is  the  amplitude  frequency  char¬ 
acteristic  of  oscillating  tube  of  force  with  allowance  for  the  resonance  properties.  (A  ss 
1-10).  The  value  of  the  height  integrated  Pedersen  conductivity  Ep  varies  from  .1  mho 
for  night  time  conditions  to  1-10  mho  during  daytime  conditions.  The  value  of  is 
of  the  order  of  25  mV/m  in  tlie  high  latitude  ionosphere  and  reaches  values  of  100-150 
inV/m  during  substorms.  It  is  obvious  that  night  time  conditions  are  favoied.  Values  of 
AEx  ~  10  mV/rn  compatible  with  threshold  require  values  of  AE  .1  mhc.  This  implies 
an  energy  deposition  into  electrons  of  10*®  eV/m^  per  ULF  p'^riod.  For  an  area  of  25  km^, 
the  required  energy  deposition  to  electrons  will  be  the  order  1-2  kJ.  Since  the  frequency 
is  of  the  order  of  a  second,  the  required  power  absorption  by  electrons  will  be  of  the  order 
of  1-2  kW.  Given  the  inefficiency  ot  long  time  electron  heating  in  the  D-E  region  it  will 
require  an  HF  facility  with  1-2  MW  ground  power  operating  in  the  4-7  MHz  region. 

C.  Key  Research  Issues 

1.  What  is  the  value  of  the  electric  field  generated  in  the  region  for  frequencies  in  the 
range  .01-10  Hz  as  a  function  of  the  HF  power  density  and  frequency  and  ambient 
conditions? 

2.  What  is  the  Q  value  (i.e.,  value  of  A)  for  the  magnetospheric  shell  cavity  modes  and 
their  harmonics? 

3.  What  is  the  conversion  efficiency  from  the  shear  to  the  kinetic  Alfven  wave  in  the 
equatorial  plane? 

4.  Are  there  any  feedback  effects  due  to  the  enhamced  precipitation? 

.4  combination  of  satellite  and/or  rocket  experiments  in  conjunction  with  the  operation 
of  a  strong  high  frequency  (i.e.,  5-12  MHz)  heater  can  provide  quantitative  answers  to 
these  questions.  Preliminary  experiments  can  be  performed  using  the  Tromso  Max  Planck 
facility.  However  the  planned  HAARP  facility  could  be  a  good  instrument  to  investigate 
the  FRIP  concept. 
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FIGURE  CAPTIONS 


1.  Time  evolution  of  various  quantities  for  u,’/Q  =  1.97,  u-’p/fl  =  0.3,  a  —  80°,  and 
e  =  0.19.  The  wave  is  an  R-x  mode.  The  unmarked  lines  in  d  and  c  are  obtained 
assuming  iq  =  nyc. 

2.  Time  evolution  of  F{6  >  10).  Time  evolution  of  various  quantities  for  uj/0,  =  1.97, 
cjp/n  =  0.3,  o  =  80°,  and  e  =  0.19.  The  unmarked  curve  is  from  the  diffusion  code 
and  the  marked  curve  is  from  the  particle  code. 

3.  Time  evolution  of  various  quantities  for  u)/Q,  =  1.97.  uipJQ  =  0.3,  a  =  80°,  and 
€  =  3.9.  The  wave  is  an  R-x  mode.  The  unmarked  lines  in  d  and  c  are  obtained 
assuming  rq  =  nyc. 

4.  Time  evolution  of  F(S  >  10).  Time  evolution  of  various  quantities  for  uj/Q,  =  1.97, 
ijjp/Q,  =  0.3,  Q  =  80°,  and  e  =  3.9.  The  unmarked  curve  is  from  the  diffusion  code 
and  the  marked  curve  is  from  the  particle  code. 

5.  Time  evolution  of  vanous  quantities  for  uj/Q  =  1.97,  u>p/Cl  =  0.3,  a  =  80°,  and 
e  =  10.0.  The  wave  is  an  R-x  mode.  The  unmarked  lines  in  d  and  c  are  obtained 
assuming  =  nyc. 

6.  Time  evolution  of  F(6  >  10).  Time  evolution  of  various  quantities  for  u;/fi  =  1.97, 
Up/Q.  =  0.3,  a  =  80°,  and  e  =  10.0.  The  unmarked  curve  is  from  the  diffusion  code 
and  the  marked  curve  is  from  the  particle  code. 

7.  D ! Dqi  as  a  function  of  wave  amplitude.  The  scale  of  the  x-axis  is  linear  but  the  scale 
from  0.19  to  0.99  is  different  than  that  fron  0.99  to  10. 

8.  Time  evolution  of  various  quantities  for  =  2.6,  Wp/fi  =  0.3,  a  =  80°,  and 

e  =  0.427.  The  wave  is  an  L-o  mode.  Some  particles  are  accelerated  in  the  negative 
z-direction  and  the  approximation  Vz  =  nyc,  given  by  the  unmarked  line,  is  no  longer 
adequate. 

9.  Time  evolution  of  various  quantities  for  ui /Q,  =  2.6,  w'p/f2  =  0.3,  n  =  80°,  and  e  =  3.9. 
The  wave  is  an  L-o  mode.  Some  particles  a  e  accelerated  in  the  negative  z-direction 
and  the  approximation  tq  =  uyc,  given  by  the  unmarked  line,  is  no  longer  adequate. 

10.  Time  evolution  of  various  quantities  for  cj/Q  =  2.6,  u/p/Q,  —  0.3,  a  =  80°,  and  e  =  10.0. 
The  wave  is  an  L-o  mode.  Some  particles  are  accelerated  in  the  negative  z-direction 
and  the  approximation  v.  =  nyc,  given  by  the  unmarked  line,  is  no  longer  adequate. 
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11.  Schematic  illustration  of  particle  motion  in  the  Earth’s  ionosphere. 

12.  The  limits  of  the  stochastic  region  of  velocity  space  for  Q  —  30.23.  The  crosses  give 
the  numerically  observed  values  [adapted  from  Ref.  23]. 

13.  Amplitude-frequency  characteristics  at  the  equatorial  plane,  assuming  a)  a  symmetric 
ionosphere,  b j  an  assymetric  ionosphere  [adapted  from  Ref.  27]. 

14.  Amplitude-frequency  characteristics  for  Ab  at  the  equatorial  plane,  assuming  a  ho¬ 
mogeneous  magnetosphere  and  a  symmetric  ionosphere  [adapted  from  Ref.  27]. 

15.  Shown  is  a  schematic  illustration  of  field  resonance  induced  precipitation.  A  standing 
kinetic  Alfven  wave  has  a  component  parallel  to  field  lines  and  resonant  with  the 
electron  bounce  motion  which  induces  precipitation. 

16.  Frequency  vs.  length  for  a  simple  field  model. 
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Nature  of  the  Diffusion  above  the  Chaotic  Threshold 


H.  Karimabadi,  C.  R.  Menyuk(“) 

Deaprtment  of  Electrical  and  Computer  Engineering 
University  of  California,  San  Diego,  California  92093 

Abstract 

The  diffusion  coefficient  for  particles  in  a  magnetized  plasma  and  in  the  presence  of  a  coher¬ 
ent  wave  was  calculated  numerically.  The  results  were  compared  with  the  solution  of  a  diffusion 
equation  based  on  the  quasilinear  theory.  The  ratio  of  the  numerical  to  the  quasilineaj  diffusion 
coefficient  DfDqi  was  found  to  be  between  0.9  and  1.4  for  wave  amplitudes  sufficiently  above  the 
stochasticity  threshold.  No  coherent  oscillation  of  DJDqi,  as  a  function  of  wave  amplitude  was 
observed.  The  above  results  held  up  in  all  of  our  parameter  study. 


(a)  Department  of  Electrical  Engineering,  University  of  Maryland,  Baltimore,  Maryland  2122S. 
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The  quasiiinear  theory^  and  its  domain  of  validity  has  remained  an  area  of  active  research  due 
to  its  relevance  to  so  many  different  areas  of  physics.  Recent  advances  in  nonlinear  dynamics  have 
revealed  that  even  systems  of  low  dimensionality  can  exhibit  comple.x  behavior.  Thus,  more  recent 
tests  of  the  quasiiinear  theory  have  been  for  systems  of  low  dimensionalities.  Using  the  standard 
map  to  model  a  chaotic  system,  several  authors-  found  that  the  true  diffusion  oscillates  about  the 
quasiiinear  value  as  a  function  of  the  nonlinearity  parameter  with  the  discrepancy  between  the  two 
reaching  as  high  as  2.5.  The  key  question  is  how  much  of  the  discrepancy  is  due  to  the  nature  of 
the  approximations  made  and  how  much  is  due  to  the  limitations  of  quasiiinear  theory.  In  this 
letter,  we  consider  the  validity  of  the  quasiiinear  theory  in  a  physical  system.  We  investigate  the 
problem  of  particle  diffusion  in  a  magnetized  plasma  and  in  the  presence  of  a  coherent  wave.  We 
derive  a  diffusion  equation  and  model  the  diffusion  coefficient  based  on  the  resonance  broadening 
theory^.  We  have  developed  an  algorithm  that  allows  a  rapid  and  accurate  solution  of  the  diffusion 
equation.  Thus  we  have  been  able  to  carry  out  a  detailed  study  of  the  particle  diffusion  and  its 
comparison  withe  quasiiinear  theory. 

In  a  magnetized  plasma,  there  are  many  resonances  given  by  u;  -  fci|Ui|  =  When  the 

wave  amplitude  is  small,  the  resonances  are  S'  parated  and  the  particle  motion  is  integrable.  .A.s 
the  wave  amplitude  becomes  larger  than  a  certain  threshold,  the  adjacent  resonances  overlap  and 
particle  motion  becomes  chaotic^  and  can  in  turn  be  desribed  by  a  diffusion  equation.  We  start  our 
analysis  with  the  usual  diffusion  equation  in  2-D  that  is  averaged  over  the  gvroangle.  The  solution 
of  the  2-D  diffusion  equation  is  in  general  difficult  and  compuf  ..tionally  e.xpensive.  It  has,  however, 
been  shown^  that  the  particles  follow  the  constant  Hamiltonian  surface  in  phase  space.  Thus  by 
transforming  to  a  coordinate  system  that  has  the  Hamiltonian  as  one  of  the  axes,  we  can  further 
reduce  the  diffusion  equation  to  1-D.  The  decorr  ion  due  to  stochasticity  is  incorporated  into 
the  diffusion  coeeficient  by  means  of  the  resonanc  oroadening  theory.  The  usual  delta  function 
associated  with  each  resonance  in  the  diffusion  coefficient  is  now  replaced  by  a  bo.x  having  a  width 
of  2>.  *;  and  a  height  proportional  to  I/u-'i,,  where  u/j  is  the  bounce  frequency.  The  resulting  diffusion 
equation  is: 


a 

dt 


i_a_ 

7^7 


07 


(la) 
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where 


In  the  regions  where  resonances  o/erlap,  the  height  of  the  boxes  adr'..  The  resulting  diff’isior. 
coefficient  is  quite  bumpy  as  shown  in  Figure  1  and  we  use  the  finite  element  technique  to  soive 
the  diffusion  equation. 

Wo  have  tested  the  code  by  running  it  for  cases  where  the  diffusion  equation  has  analytic 
solutions  and  found  the  agreement  between  the  results  of  simulation  and  analysis  to  be  excenent. 
The  number  of  particles  are  typically  constant  to  better  than  one  part  in  lO-®. 

In  order  to  test  the  quasilinear  theory,  we  have  solved  the  diffusion  equation  for  many  different 
parameters  and  compaxed  the  results  with  exact  orbit  calculations.  We  take  the  initial  particle 
distribution  function  to  be  a  delta  function  centered  at  7  =  1.  In  the  exact  orbit  calculations,  we 
follow  the  -rbits  of  256  pcirticles,  all  with  initial  7=1  and  phaises  distributed  uniformly  between  0 
and  2~,  up  to  fit  =  3000.  The  reason  for  the  choice  of  the  integration  time  is  made  clear  below.  The 
resulting  oustribution  function  is  very  bumpy  due  to  the  small  number  of  particles  used.  To  obtain 
smooth  distributions,  millions  of  particles  would  need  to  be  used  -  an  approach  which  is  both  both 
impractical  and  unnecessary.  Instead,  we  compare  the  moments  of  the  two  distribution  functions. 
The  small  number  of  particles  which  we  use  is  sufficient  In  this  case  to  yield  an  accurate  result.  In 
what  follows,  we  use  the  quantities  <7—70  >=<  A7  >  and  <  (7  —  7o)'  >  =  <  >  to  make 

the  comparison,  where  70  is  the  initial  7  of  the  particles  and  the  brackets  indicate  an  average  over 
the  initial  conditions: 

<  (7  -  lof  >  = 

Note  that  Dy-j  a<  (^7)^  >. 

We  present  the  results  for  two  sets  of  parameters:  First,  we  consider  the  diffusion  due  to  a 
R-x  mode  with  w/fl  =  1.97,  a  =  80°,  and  =  0.3.  The  values  of  stochasticity  threshold 

Cihrs  calculated  from  analysis  and  simulations  are  0.19  and  0.143,  respectively.  We  carry  out  our 
comparison  for  wave  amplitudes  ranging  from  0.19  to  10.  The  variation  of  several  quantities  as 
functions  of  time  are  snown  in  Fig.2  for  e  =  0.19.  The  quantity  <  (^7)"  >  is  seen  to  follow  a  linear 
increase  in  time  in  accord  with  the  underlying  assumptions  built  into  the  Fokker-Planck  equation. 
The  particles  have  not  had  enough  time  to  sample  all  the  available  resonances  and  as  a  result  the 


/”(7-7o)V(7.0<^7 

/r/(7,0<^7 
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quaniity  <  (dy)'  >  increases  with  time.  Once  the  particles  reach  the  energies  beyond  which  D-,. 
becomes  very  sraaii.  the  <  (^7)'  >  flattens  out  as  a  function  of  time  and  diffusion  stops.  For  the 
parameters  of  Fig. 2,  the  flattening  out  does  not  happen  until  Of  3000. 

The  quajitity  D f  Dql  which  is  the  ratio  of  the  e.xact  orbit  calculation  of  <  (Aj)'  >  to  that 
obtained  from  the  diffusion  code,  deviates  significantly  from  unity  at  early  times  (Fig. 10c)  but  settles 
down  close  to  its  asymptotic  value  for  Clt  ^  750-  This  is  expected  since  the  diffusion  formalism  is 
strictly  valid  for  times  long  compared  with  correlation  time  scales.  Similar  results  were  also  seen 
at  higher  amplitudes. 

The  deviation  of  particle  motion  from  the  predicted  quasillnear  diffusion  at  early  times  can 
occur  due  to  at  least  two  effects:  First,  there  can  exist  small  but  finite  islands  of  stability  within 
an  otherwise  stochastic  phase  space.  A  particle  coming  close  to  such  an  island  can  get  temporarily 
trapped.  Such  stickiness  in  phase  space  can  obviously  have  an  effect  on  the  diffusion  of  the  particles, 
.ks  it  turns  out,  however,  this  effect  is  not  that  important  in  general  and  does  not  play  a  significant 
role  for  9.t  >  750.  As  we  mentioned  earlier,  the  trapping  width  increases  as  a  function  of  harmonic 
number  and  in  the  transition  to  global  chaos,  the  border  of  chaos  lies  at  the  lowest  energies.  Since 
we  start  the  particles  at  7  =  1,  and  the  stickiness  is  most  important  at  the  border  of  chaos  which 
is  also  at  7  ~  1,  the  particles  often  experience  some  initial  stickiness;  however,  once  the  particles 
reach  the  higher  energies,  they  diffuse  freely.  Thus,  there  exists  a  finite  time  before  all  the  particles 
can  diffuse  freely.  The  second  source  of  deviation  from  the  predicted  quasillnear  behavior  is  that 
the  particles  may  initially  sample  a  few  resonances  in  which  case  the  motion  is  closer  to  a  coherent 
acceleration  than  diffusion.  The  first  effect  results  in  a  retardation  in  the  diffusion  and  the  second 
effect  leads  to  diffusion  larger  than  Dqi-  The  balance  between  these  two  effects  determines  £)/ Dqi 
at  early  times.  For  Q.t  >  750,  the  above  two  effects  become  much  less  important. 

We  have  plotted  <  7  >  /  <  7^  >  and  DfDqi  at  9.t  =  3000,  for  a  range  of  wave  amplitudes 
as  shown  in  Fig. 3.  The  quasillnear  theory  appears  to  be  quite  adequate  in  describing  the  particle 
motion  over  this  rather  large  range  in  wave  amplitude.  For  e  ^  10,  the  particle  motion  approaches 
the  integrable  regime  of  unmagnetized  plasma  in  which  case  the  diffusion  equation  would  no  longer 
be  valid.  The  deviation  of  DjDqi,  from  unity  as  a  function  of  wave  amplitude  has  a  sporadic 
behavior  and  does  not  seem  to  follow  a  nice  oscillating  pattern.  In  short,  we  have  found  no  evidence 
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for  oscillations  of  Dj Dqi  as  a  function  of  wave  amplitude,  in  contrast  to  earlier  studies  based  on 
the  standard  mapping*.  Furthermore,  the  fluctuation  amplitude  of  DjDqc  is  much  smaller  than 
that  seen  in  previous  numerical  work.  The  presence  of  resonances  of  different  sizes  ciearly  leads  to 
averaging  beyond  that  which  is  found  in  the  standard  map  and  thus  to  a  closer  correspondence  to 
quasilinear  theory.  These  results  suggest  that  the  standard  map  is  not  an  ideal  paradigm  for  real, 
physical  systems. 

Next,  we  consider  the  diffusion  eqn  due  to  the  L-o  mode  with  w/fl  =  2.6,  q  =  80°,  and 
utp/fl  =  0.3.  The  value  of  threshold  based  on  analysis  is  0.42,  and  the  wave  amplitude  necessary  to 
accelerate  1%  and  60%  of  the  particles  stochastically  from  zero  energy  is  0.4  and  0.75,  respectively. 
The  quantity  <  A7  >  /  <  A7  >qi  is  well  below  unity  for  c  =  0.43  and  it  does  not  reach  unity 
until  the  wave  amplitude  has  been  increased  to  t  ~  1.5.  For  c  ^  1.5,  the  <  A7  >  /  <  A7  >qL  has 
small  fluctuations  about  unity.  The  <  A7  >ql  was  obtained  by  using  H  =  1.  Phase  effects  are 
important  for  the  L-o  mode  and  the  Cthn  depends  strongly  on  the  value  of  H  which  has  a  finite 
range  even  for  7=1.  The  Cthr*  can  change  by  more  than  a  factor  of  3  depending  on  the  value  of  //. 
Thus,  we  must  take  into  account  the  fraction  of  particles  that  are  not  stochastic  at  a  given  wave 
amplitude.  The  ratio  of  stochastic  particles  is  ~  60%  at  e  =  0.43,  ~  75%  at  e  =  0.8,  ~  80%  at 
«  =  1.1,  ~  93%  at  f  =  1.3,  and  ~  99%  at  e  =  1.5.  This  is  why  <  A7  >  /  <  A7  >Qi,  is  below  unity 
for  e  1.5.  Let  us  denote  the  integrable  particles  by  the  subscript  ’int’  and  the  stochastic  particles 
by  the  subscript  ’st’.  What  we  should  really  be  plotting  in  Fig. 4a  is  <  A7  >  ,c  /  <  A7  >qL  and 
not  <  A7  >  /  <  A7  >qL.  But  <  7  >*<  /  <  7  >QL~  <  1  >  I  <  -(  >qL,  where  we 

have  assumed  <  7  >int  /  <  7  ><?£,<  1-  Using  the  values  of  N,t  cited  above  for  various  e’s,  it  then 
follows  that  <  A7  >  ,t  I  <  A7  >qL  is  indeed  very  close  to  unity.  In  short,  the  discrepancy  between 
<  A7  >  and  <  A7  >Qr,  in  Fig.4a  for  c  ;$  1.5,  is  due  to  the  fact  that  we  have  included  integrable 
orbits  in  <  7  >  but  not  in  <  A7  >qL-  The  fraction  of  stochastic  particles  can  be  easily  calculated 
analytically  by  incorporating  phase  effects  in  the  expression  for  the  stochasticity  threshold.  W'e 
then  find  that  stochastic  particles  are  well  described  by  the  diffusion  code. 

The  agreement  between  D  and  Dq^  is  again  very  good  for  a  large  range  of  wave  amplitudes 
as  shown  in  Fig.4b.  This  is  encouraging  considering  that  the  theoretical  analysis  presented  in  this 
paper  are  at  their  worst  for  the  parameters  in  Fig.4  where  the  phase  effects  are  important  and  H 
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i5  no:  just  confined  to  one  value.  Thus,  we  conclude  that  the  diffusion  formalism  is  quite  robust 
and  is  highly  accurate  in  predicting  the  time  evolution  of  an  ensemble  of  particles. 
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Figure  Captions 


Figure  1.  Diffusion  Coefficient.  The  parameters  used  are  Wp/D  =  0.3,  w/fi  =  1.97,  a  =  S0°, 
and  e  =  0.19.  The  wive  is  a  R-x  mode.  Only  the  diffusion  coefficient  between 
7  =  1  and  7  =  40  is  shown. 

Figure  2  Time  evolution  of  various  quantities  for  a;/fi  =  1.97,  Wp/0  =  0.3,  a  =  80°,  and 
€  =  0.19.  The  wave  is  a  R-x  mode. 

Figure  3.  The  variation  of  <  A7  >  /  <  A7  >q£,  and  D/Dql  as  functions  of  wave  amplitude. 

The  waves  are  R-x  modes  with  tj/fi  =  1.97,  a  =  80°,  and  Wp/D  =  0.3.  The  scale  in 
the  x-axis  is  linear  but  the  scale  from  0.19  to  0.99  being  different  than  that  from 
0.99  to  10. 

Figure  4.  The  variation  of  <  A7  >  /  <  A7  >(}£,  and  DIDql  as  functions  of  wave  ampli¬ 
tude.  The  waves  are  L-o  mode  with  uj/Q,  =  2.6,  a  =  80°,  and  uip/Cl  =  0.3.  The 
discrepancy  between  <  A7  >  and  <  A7  >ql  for  e  ^  1.5  is  due  to  the  presence 
of  integrable  orbits.  If  we  exclude  the  integrable  orbits  in  calculating  <  A7  > , 
we  once  again  find  an  excellent  agreement  with  the  diffusion  code.  The  fraction  of 
integrable  orbits  can  be  determined  apriori  analytically.  The  scale  in  the  i-a.xis 
is  linear  but  the  scale  from  0.43  tp  1.5  is  different  than  that  from  1.5  tp  10. 
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Abstract 

Electron  acceleration  by  a  wave  propagating  obliquely  in  a  magnetized  plasma  is  considered. 
For  wave  amplitudes  above  the  stochasticity  threshold,  the  motion  of  electrons  is  diffusive.  A  test 
of  a  simple  diffusion  formalism  is  presented.  We  reduce  the  diffusion  equation  to  a  1-D  form  by 
transforming  to  a  coordinate  system  that  has  the  Hamiltonian  as  one  of  the  a.xes.  We  present  a 
model  for  the  diffusion  coefficient  and  solve  the  resulting  diffusion  equation  by  means  of  the  finite 
element  technique.  The  results  are  compared  with  numerical  solutions  of  the  orbits.  Finally,  we 
apply  our  results  to  the  problem  of  electron  diffusion  in  the  lower  ionosphere. 

(a)  Permanent  Address:  Department  of  Electrical  Engineering,  University  of  Maryland,  Baltimore, 
Maryland  21228. 
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I.  Introduction 


The  problem  of  charged  particle  acceleration  by  waves  is  of  fundamental  importance  in  plasma 
physics  and  plays  a  central  role  in  the  understanding  of  many  processes  in  space  physics.  In  an 
unmagnetized  plasma,  the  motion  of  a  charged  particle  interacting  with  a  plane  wave  is  very  simple: 
there  exist®  at  most  one  resonance  and  the  system  is  integrable.  The  trapping  width  is  proportional 
to  where  t  is  the  wave  amplitude.  The  particle  motion  can  then  become  random  only  if  there 
exist  several  waves  with  amplitudes  and  frequencies  such  that  the  trapping  width  of  adjacent  waves 
overlap.  In  other  words,  the  condition  for  the  particles  to  be  stochastic  is  for  the  bounce  frequency 
in  the  potential  well  of  each  wave  to  be  larger  than  the  spacing  between  the  adjacent  resonances 
Ub  >  Aw.  The  evolution  of  an  ensemble  of  particles  can  then  be  described  by  the  quasilinear 
(diffusion)  equation  [Vedenov,  Velikhov,  and  Sagdeev,  1961;  1962;  Zaslavsky,  and  Chirikov,  1972]. 
Conversely,  when  the  wave  amplitudes  become  large  enough  to  satisfy  the  condition  wj  >  N^Au, 
then  the  stochasticity  disappears  as  the  effect  of  individual  resonances  would  be  washed  out  and 
the  particle  would  see  a  single  potential  weU  with  slowly  varying  parameters.  Here  N ^  is  the 
number  of  waves  in  the  wave  packet.  Thus,  the  condition  for  the  validity  of  the  quasilinear  theory 
is  N^Au  >  Wf,  >  Aw. 

The  presence  of  a  static  magnetic  field  changes  the  particle  motion  in  two  important  ways; 
(1)  particles  can  resonate  with  the  harmonics  of  the  gyrational  motion,  (2)  resonances  can  occur 
even  if  the  wave  has  a  phase  velocity  larger  than  the  speed  of  light.  The  resonance  condition  in  a 
magnetized  plasma  is  w  -  /:||U||  =  where  f  =  0,  ±1,  ±2, ...,  and  fl  =  \q\BolTnc.  The  trapping 

width  associated  with  each  resonance  is  proportional  to  As  long  as  c  1  the  neighboring 

resonances  are  well  separated  and  the  particle  motion  is  periodic.  At  a  critical  value  of  the  wave 
amplitude  referred  to  as  the  stochasticity  threshold  [e.g.,  Lichtenberg  and  Lieberman,  1983],  the 
adjacent  resonances  overlap  and  the  motion  becomes  random  or  chaotic.  The  particle  can  then 
sample  several  resonances  and  gain  large  energies  in  the  process.  In  this  case,  the  motion  of  the 
particle  can  still  be  described  by  the  quasilinear  theory  even  though  there  is  only  one  wave  present. 

The  general  problem  of  single  particle  raotiori  under  the  influence  of  an  obliquely  propagating 
plane  wave  has  been  treated  in  a  series  of  papers  by  Karimabadi  and  collaborators  [Menyuk.  el 
al.,  1987;  1988;  .Akimoto  and  Karimabadi,  1989;  Karimabadi,  et  al.,  1990].  The  main  result  of 
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these  single  particle  studies  has  been  the  realization  that  for  waves  with  .Vi|  <  1.  the  Hamiltonian 
surfaces  become  topologice_iy  open,  thus  allowing  large  particle  acceleration.  Here.  A'.i  is  the 
refractive  index  parallel  to  the  magnetic  field.  In  this  paper,  we  derive  a  diffusion  equation  that 
describes  the  evolution  of  particle  distribution  function  for  wave  amplitudes  above  the  stochasticity 
threshold  and  A’jj  <  1.  The  generalization  of  our  work  to  the  case  of  Ahj  >  1  is  trivial.  We  reduce 
the  diffusion  equation  to  1-D  by  transforming  to  a  coordinate  system  with  the  Hamiltonian  as  one 
of  the  axes.  We  then  present  a  model  for  the  diffusion  coefficient  and  solve  the  diffusion  equation 
using  the  finite  element  technique.  This  algorithm  is  found  to  be  very  fast  and  accurate.  The 
solutions  of  the  diffusion  equation  are  compared  with  exact  numerical  solutions  of  the  orbits.  The 
conditions  under  which  the  diffusion  formalism  breaks  down  are  Discussed. 

The  techniques  and  algorithm  e.xpounded  in  this  paper  are  quite  general  and  are  likely  to  find 
applications  in  many  problems  in  space  physics  and  astrophysics.  We  have  already  applied  the  above 
techniques  successfully  to  the  problem  of  par.icle  acceleration  in  cometary  shocks  [ivarimabadi,  et 
ai,  1989]  and  in  supernova  remnants  [Kaiimabadi,  1988).  In  section  [VI],  we  briefly  apply  our  results 
to  the  problem  of  remote  acceleration  of  ionospheric  electrons  from  ground  based  radio  transmitters. 
Accelerated  electrons  could  be  used  to  create  an  artificial  aurora,  to  probe  the  magnetosphere,  and 
to  measure  the  cross  section  of  ionospheric  constituents  interacting  with  energetic  electrons.  We 
find  the  power  needed  for  stochastic  electron  acceleration  to  be  beyond  the  present  day  capabilities. 
We  propose  ways  to  overcome  this  problem.  We  present  our  conclusion  in  section  fVIIl. 

II.  Test  Particle  Results 

In  this  section  we  describe  the  formulation  of  the  resonance  overlap  criterion  which  w-e  use  to 
determine  the  acceleration  threshold.  More  details  can  be  found  in  Karimabadi.  et  al.  [19901. 

In  general,  we  assume  that  each  electron  is  moving  in  a  homogeneous  magnetic  field  and  a  plane 
wave  propagating  obliquely  to  the  magnetic  field.  This  plane  wave  can  have  both  electromagnetic 
and  electrostatic  components.  We  thus  write  the  particle  Hamiltonian  as 


H  = 


[cP-qA)' 


+  m*c‘* 


1/2 


+  q^  —  mc~~i  +  . 


la) 


where 


A  =  .4i  ~  sin  cm-  +  Aj  cos  ihcy  -  .4i  sin  gm.  +  xD.,Cy  =  A,m  m  xB^ey,  (lb) 


$  =  <J>o  sin 


(Ic) 


and 


Tp  =  ki_x  +  ^(1^  —  ijJt. 


(Id) 


Although  we  aje  primarily  interested  in  electron  acceleration,  we  keep  the  charge  9  arbitrary  in  the 
derivation  of  the  equations.  We  next  use  a  series  of  canonical  transformations  and  the  techniques 
of  Hamiltonian  perturbation  theory  to  reduce  the  Hamiltonian  to  the  form 


(2) 


The  zero-order  Hamiltonian  is  given  by 


A  LJ 

Ho  =  mc^jo  -  T-Pii, 


'3a) 


where 


70  =  (  1  +  ^2  2 


_Plll 

rn^c^ 


1/2 


(3b) 


The  quantities  px  and  py  represent  the  momentum  perpendicular  and  paraJlel  to  the  magnetic 
field.  The  surfaces  of  constant  Ho  are  elliptic  in  p||  -  px  space  if  Afu  >  1,  hyperbolic  if  fV||  <  1  and 
parabolic  if  jV||  =  1 .  Equation  (2)  is  derived  by  expanding  in  powers  of  the  wave  field;  it  is  consistent 
with  this  expansion  to  takepx  and  py  to  be  purely  mechanical.  The  first-order  Hamiltonian  is  given 
by 

00 

Hi  -  ^  Z(S\vi{k\\z  +  19)  ,  (4) 
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/here  9  is  the  gyroangle  and 


me"  I 
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Q  ^  A 

— - sin  a  -  —  - - f.  cos  a 

7  j//|  fcxc 
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Cl  = 
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\(i\A2  \q\E2 


=  =  ,5d) 

mc^  mcNu 

Here,  the  quantity  a  is  the  angle  of  wave  propagation  with  respect  to  the  magnetic  field,  q  is  the 
electric  charge,  and  fl  is  the  magnitude  of  the  nonrelativistic  gyrofrequency.  The  explicit  time 
dependence  in  the  original  Hamiltonian  of  Eq.  (la)  has  been  absorbed  into  the  2 -coordinate  of 
Eq.  (4). 

Resonances  occur  when  A:||i  +  =  0,  or 

=  (6) 

rr^lo  Jo 

a  result  which  can  also  be  obtained  from  standard,  linear  theory.  It  can  be  easily  verified  that  the 
resonance  surfaces  are  elliptic  (hyperbolic)  when  the  /T-surface  is  hyperbolic  (elliptic).  Typically, 
for  a  given  value  of  C  and  H,  there  will  be  only  one  pair  (px,P(|)  which  satisfies  Eq.  (6).  It  is  this 
case  which  we  are  interested  in  treating.  In  some  cases,  it  is  possible  for  Eq.  (6)  to  be  satisfied  for 
every  point  (pxiP||)  on  the  Hamiltonian  surface  in  phaise  space.  In  these  cases,  special  methods  of 
treatment  are  required  [Karimabadi  et  al.,  199UJ.  We  will  not  discuss  this  issue  in  this  paper. 

If  the  wave  amplitude  is  smsdl,  then  an  electron  will  be  affected  by  at  most  one  resonance  In 
this  case,  all  but  that  single  resonance’s  contribution  to  the  Hamiltonian  can  be  ignored,  and  the 


Hamiltonian  becomes 


//  =  Ho(p\\,Pi  )  +  Z,sin(fc||2), 


where  z  has  been  redefined  to  absorb  id.  In  this  re-defined  coordinate  system,  the  resonance  occurs 
at  a  point  fp||^,pxr)  where  i  =  0.  Writing 


letting 


1  0^  H 

-  ^o{P\\r.PXr)  +  2~^ 


(T’II  -  /hlr) 
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and  absorbing  the  constant  HJp]i^r,Pir)  into  H,  the  Hamiltonian  becomes 


10) 


which  is  the  pendulum  Hamiltonian.  The  trapping  half-width  and  the  bounce  frequency  are  then 
given  by 

Apii  =  (11) 


and 


,1/2 


M 


CJf,  =  ^11 

respectively.  The  trapping  width  in  7  can  be  easily  found  from  Apy: 

^Pll 


A7  = 


mcNii 


(12) 


(13) 


The  trapping  half-width  given  by  (11)  is  independent  of  the  initial  phase.  This  is  due  to  our 
approximation  of  the  Hamiltonian  near  a  resonance  by  the  pendulum  Hamiltonian.  The  phase 
dependence  of  the  trapping  half-width  is  mainly  due  to  the  py  component  of  momentum  and  the 
electrostatic  component  of  the  wave.  If  c  0,  the  phase  dependence  of  the  trapping  half-width 
becomes  increasingly  important  as  q  — ►  90°.  The  above  anailysis  can  be  easily  generalized  to  include 
the  effect  of  phase  on  the  trapping  width  [Ginet  and  Heinemann,  1990].  One  then  ends  up  with 
an  equation  for  the  trapping  half-width  that  can  in  general  only  be  sol/ed  numerically.  This  is  not 
very  useful  to  us  and  w’e  use  (11)  in  this  paper. 

It  is  important  to  note  that  the  term  “trapped”  has  a  very  specific  meaning:  it  refers  to  orbits 
close  enough  to  resonance  condition  so  that  the  sum  over  all  the  harmonics  can  be  replaced  by  a 
single  term  in  the  sum.  The  particles  not  very  close  to  the  resonance  can  still  be  accelerated,  but 
the  acceleration  is  much  less  effective. 

In  deriving  (11)  and  (12),  we  have  explicitly  assumed  that  a  ^  90°.  For  q  =  90°,  most  of 
the  steps  that  led  to  (11)  remain  the  same.  The  only  difference  comes  in  eliminating  the  time 
dependence  in  the  Hamiltonian:  the  transformation  F2  =  P^||[‘  —  wd/A:||]  is  replaced  by  = 
p\( -.jt/k).  It  is  then  easy  to  show  that  the  trapping  width  is  still  given  by  (13).  The  only 
difference  is  that  for  q  =  90°,  py  is  a  constant  of  motion  and  is  determined  a  priori  from  the  initial 
conditions.  The  sensitivity  of  the  trapping  width  to  py  and  thus  to  the  initial  phase  at  a  ~  90°  is 
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not  surprising  since  the  term  £(sin  v)/j'''’||  in  H  becomes  very  iajge  as  cy  —  90°.  The  Hamiltonian 
of  a  particle  initially  at  rest  can  then  have  a  large  range  of  values  depending  on  the  initial  phase. 
For  Q  ~  90°,  p||  is  still  given  by  the  solution  of  (3a)  and  (6)  but  with  Ho  replaced  by  H  ana  Vv. 
replaced  by  the  canonical  momentum.  If  fi  =  0,  the  trapping  width  becomes  less  sensitive  to  the 
initial  phase. 

The  septiration  on  an  /f— surface  between  two  adjacent  primary  resonances  is 

mc^Q.k\\ 


<Jpil  = 


The  stochasticity  threshold  condition  based  on  the  overlap  criterion  [Chirikov.  1979]  then  reads 


14) 


AP|| 

=  4  a;j/(n/7o)  ^  1- 


(15) 


As  long  as  the  resonances  do  not  overlap,  an  electron  in  the  neighborhood  of  the  resonance  i  can 
change  in  p||  by  at  most  2Ap||.  Once  resonances  overlap,  the  electron  can  jump  from  resonance 
to  resonance  and  in  this  way  gain  large  amounts  of  energy.  VVe  can  determine  in  this  way  the 
threshold  wave  intensity  to  accelerate  electrons  originally  at  rest  to  10  Mev  or  more. 


III.  The  Diffusion  Equation 

There  are  two  equivalent  ways  of  deriving  the  diffusion  equation.  The  first,  where  the  ef¬ 
fects  of  particle  motion  on  the  imposed  fields  are  neglected,  lends  tc  the  Fokker-Planck  equation 
[Chandrasekhar,  1943;  Romanov  and  Filippov.  ICGl;  Sturrock,  1966] 


tl6a) 


where 


D  =  lim  — ; — , 
*00  At 


!  16bl 


D  = 


lim 

*00 


(UW) 
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!  1 6c  I 


and  where  the  brackets  indicate  an  average  over  initial  conditions.  In  the  derivation  of  I  16a),  it  is 
assumed  that:  (1)  The  change  in  velocity  Au  in  a  ti...c  interval  At  is  such  that  (Au)  and  (Ac‘) 
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contain  contributions  linear  in  At:  whereas,  all  higher  products  have  expectation  values  that  vary 
as  a  higher  power  of  At  and  are  therefore  neglected.  (2)  The  trajectory  which  a  particle  follows 
depends  only  on  the  instantaneous  values  of  its  physical  paxameters  and  is  entirely  independent  of 
its  whole  previous  history  (Markovian  assumption). 

The  second  derivation  of  the  diffusion  equation  is  based  on  the  quasilinear  theory  [Vedenov, 
Velikhov,  and  Sagdeev,  1961,  1962;  Drummond,  and  Pines,  1962],  They  conclude  that: 


dt  dv^  dv’ 


(17) 


The  diffusion  equations  (16a)  and  (17)  are  equivalent  since  the  friction  and  diffusion  coefficients  in 
the  Fokker- Planck  equations  are  related  by 


B  = 


2  dv  ' 


(18) 


The  Hamiltonian  (la)  is  a  function  of  py  and  px  and  it  is  thus  convenient  to  start  with  the  relativistic 
diffusion  equation  in  cylindrical  coordinates 


where 
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=  Dxii  =  /  dr(px(0P||(«  +  r)). 

Jo 

In  120),  the  angular  brackets  indicate  an  average  over  the  initial  coordinates: 

Jo 


(20c) 


d9o 


:2ii 


2?r  Jo 

Equation  (19)  is  averaged  over  the  cylindrical  angle  <t>  which  becomes  randomized  on  timescales 
much  faster  than  the  relaxation  times  of  distributions  in  pu  and  px- 
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Integrating  over  the  unperturbed  orbits  we  obtain 

t  f: 

"  t=-oo  ^  11/  l=_oo 

l=-oo^^  ^  '  '  l-~oo 

Dl\\  =  £>||i  =  “  ~)  =  £ 

e=-oo^-^  ^  ^  ^  '  r=-oo 

where  S(x)  is  the  usual  Dirac  ^-function.  These  diffusion  coefficients  are  not  useful  as  they  stand 
because  in  their  derivation  we  did  not  take  into  account  the  decorrelation  due  to  stochasticity.  In 
order  to  overcome  this  problem,  we  use  a  modified  form  of  resonance  broadening  theory  [Dupree, 
1966]  as  given  below. 

A.  Diffusion  Model 

The  d-function  in  the  diffusion  coefficients  came  about  when  we  substituted  unperturbed 
orbits  in  the  dght  hand  side  of  (20),  which  then  implies  incorrectly  that  if  the  particle  is  in  resonance 
with  the  wave,  it  will  maintain  resonance  indefinitely.  We  showed  in  previous  sections  that  the 
resonance  width  is  indeed  finite  and  is  given  by  (11).  The  finiteness  of  the  trapping  width  is  due  to 
the  fact  that  A;||t;||  ^  EQ-y/-y'^,  and,  thus,  the  particle  is  shifted  from  exact  resonance  by  the  wave. 

It  can  be  easily  shown  that  the  decorrelation  time  or  the  shift  in  the  resonance  condition  due 
to  the  finite  trapping  width  is  roughly  2lji,  in  the  f-th  resonance.  Thus,  we  replace  the  6 -function 


by  a  function  /(z),  where 


_  j  1/4^6, 
10^ 


if  |i|  <  2w(, 

|i|  >  2ui, 


The  quantity  uji,  corresponds  to  the  bounce  frequency  (12)  in  the  ^th  resonance.  Note  that 


/  +  00  r  +  co 

f{x)dx  =  /  6{x)dx  =  1. 

-  OO  —  'X) 


This  spread  is  equivalent  through  the  order  to  which  we  are  working  to  the  spread  in  p  ,]  calculated 
in  (11). 


B.  Reduction  to  1-D 

Because  H  is  a  conserved  quantity,  the  diffusion  proceeds  along  a  constant  //-surface.  It 
is  then  useful  to  replace  the  (p||,px)  coordinate  system  with  the  coordinate  system,  which 
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reduces  the  diffusion  equation  to  one  dimension.  The  details  are  given  in  the  appendix  and  the 
diffusion  equation  reduces  to 


dt 


iA 

7  dj 


(23a) 


where 


The  sum  is  over  regions  of  stochasticity.  In  other  words,  we  first  calculate  a  diffusion  coefficient 
for  each  harmonic,  and,  then,  in  regions  where  resonances  overlap,  we  sum  the  contribution  due  to 
each  resonance.  We  stress  that  the  problem  is  still  two-dimensional.  We  have  reduced  the  problem 
to  one  dimension  along  each  H  —surface,  but  we  still  have  to  consider  solutions  for  many  values 
if  the  initial  distribution  function  spans  more  than  one  value  of  H. 

The  way  in  which  we  have  defined  implies  that  its  derivative  is  zero  everywhere  except 
on  a  countable  set  of  points  where  it  becomes  infinite.  Hence,  the  numerical  evaluation  of  (23) 
requires  great  care.  We  deal  with  this  difficulty  by  using  a  finite  element  discretization  method  in 
which  case  the  problem  should  be  characterized  by  a  variational  formulation.  Below,  we  use  the 
Ritz-Galerkin  method  [Muschietti,  Appert  and  Vaclarik,  1982;  Dhatti  sind  Touzot,  1984]  to  give 
the  diffusion  equation  an  equivalent  variational  formulation. 


C.  Discretization  of  the  Diffusion  Code 

We  start  by  multiplying  (23a)  by  7  and  a  wighting  function  ^(7)  and  then  integrate  over  the 
energy  domain  T  =  (7min>7max)  to  obtain 
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Of  (  df^ 


’^min  ^ 


n 


(24) 


Particle  conservation  implies  that  the  first  term  on  the  right  hand  side  of  (24)  is  zero.  Dividing  the 
range  (7min7  7max)  into  N  equal  sized  elements,  the  distribution  function  may  be  written 


25) 


j=i 


where  the  1:7(7)  are  the  basis  functions  for  the  finite  elements,  defined  so  that 


ip,  = 


(7  -  7.-i)/(7.  -  7.-i).  7.-1  <  7  <  7. 
(7i4-i  -  7)/(7.  +  i  -  7.),  7.  <  7  <  7.-t-i. 
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and  is  zero  elsewhere,  where  (7,_i,7,)  is  the  :-th  element.  We  now  use  the  v,  os  the  weishtir.z 
elements  in  (25)  and  defining  the  elements  of  the  matrices  A  and  B 


where  i  and  j  range  from  1  -  jV,  we  obtain  the  equation 


A- 


dt 


(26a) 


(26b) 


(27) 


Time  is  discretized  using  a  leapfrog  scheme,  and  (27)  is  easily  solved  as  both  A  and  B  are  tridiagonal. 
Note  that  the  matrices  A  and  B  are  time  independent  as  the  fields  are  constant  in  time  and  thus 
these  matrices  need  to  be  evaluated  just  once.  The  boundary  conditions  are  that  =  0  at  the 
lower  boundary  and  /  =  0  at  the  upper  boundary.  We  choose  7mai  to  be  sufficiently  large  that  no 
diffusion  to  this  boundary  occurs  over  the  time  scale  of  the  simulation. 

The  diffusion  coefficient  D^-,  as  given  by  (23b)  is  evaluated  numerically  as  a  function  of  as 
follows:  First,  we  determine  all  the  resonances  that  intersect  the  //-surface  in  the  interval  from 
Tmin  to  7ma*  and  Calculate  the  trapping  width  for  each  resonance.  We  model  the  contribution  of 
each  resonance  to  the  diffusion  coefficient  as  a  box  having  a  width  equal  to  the  trapping  width  and 
a  height  given  by  the  quantity  inside  the  summation  sign  in  (23b).  So,  there  is  a  box  centered 
around  each  resonance  with  known  height  and  width.  Next,  we  add  the  heights  in  regions  where 
resonances  overlap,  otherwise  leaving  the  boxes  intact.  The  resulting  diffusion  coefficient  is  not 
smooth  (Fig.  1)  and  is  not  differentiable.  However,  as  evident  from  (24)  and  (26b),  what  is  needed 
is  the  integral  of  over  each  finite  element  which  is  well  defined  and  easily  obtained. 

The  stability  condition  for  our  code  is  where  h  is  the  grid  spacing.  In  other 

words,  the  time  step  used  in  solving  the  diffusion  equation  should  be  smaller  than  the  lime  it  takes 
the  particle  to  diffuse  across  the  grid. 

We  have  tested  the  code  by  running  it  for  cases  where  the  diffusion  t-ipiation  has  analytic 
solutions  and  found  the  agreement  between  numerical  and  analytical  results  to  be  excellent.  We 
have  also  found  the  number  of  particles  to  be  typically  constant  to  better  than  one  part  in  10 '7 
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Before  we  compare  the  solutions  of  the  diffusion  code  with  exact  or'oit  calculations,  it  is  im¬ 
portant  to  check  the  accuracy  of  our  analytical  expressions  for  the  trapping  width  and  stochasticiiy 
threshold.  This  is  done  below. 


IV.  Threshold  Calculation 

Here  we  are  interested  in  determining  the  conditions  under  which  the  electrons  in  the  ionosphere 
can  be  accelerated  to  large  energies  from  rest.  The  analytical  threshold  is  based  on  the  overlap 
condition,  and  there  is  no  guarantee  that  at  the  overlap  amplitude  the  resonances  extend  down  to 
7  =  1.  Thus,  the  quantity  we  are  interested  in  is  not  the  amplitude  at  which  the  first  two  resonances 
overlap  but  the  amplitude  at  which  the  first  three  resonances  overlap  with  the  stochastic  region  in 
phase  space  reaching  do./n  to  7  =  1.  We  require  the  overlap  of  the  first  three  resonances  because 
the  waves  we  are  considering  have  rV  ~  1  and  for  such  waves  the  trapping  width  is  an  increasing 
function  of  the  harmonic  number  up  to  some  critical  t  beyond  which  the  trapping  width  drops  off 
exponentially.  We  have  found  empirically  that  if  the  first  three  resonances  overlap,  the  remaining 
resonances  will  overlap  as  well.  Thus,  the  overlap  of  the  first  three  resonances  with  the  additional 
constraint  that  trapping  width(s)  reach  down  to  7  =  1  ensures  global  stochasticity  in  the  phase 
space  region  of  interest. 

The  experimental  value  of  wave  amplitudes  for  which  zero  energy  particles  can  be  stochastically 
accelerated  are  determined  by  solving  the  orbits  of  256  particles,  all  with  initial  7  =  1  and  initial 
phases  uniformly  distributed  between  0  and  2~.  The  orbits  are  typically  solved  up  to  =  1000 
and  the  maximum  energy  7ma»  reached  for  each  particle  during  this  time  interval  is  plotted  as  a 
function  of  the  initial  phase,  as  shown  in  Fig.  2.  This  method  yields  a  quantitative  measure  of 
the  fraction  of  particles  that  are  stochastic.  In  other  contexts,  the  surface  of  section  technique  is 
useful  in  determining  the  level  of  chaos  in  phase  space;  however,  since  each  phase  corresponds  to  a 
different  value  of  //,  many  surface  of  section  plots  are  needed  to  determine  the  phase  dependence 
and,  thus,  this  is  not  practical  for  us. 

The  comparison  of  the  analytical  and  numerical  conditions  are  given  below  for  various  plasma 
parameters  and  wave  modes. 
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(i)  R-x  mode,  Wp/fi  =  0.3 

First,  we  consider  the  R-x  mode  in  a  plasma  with  Wp/fi  =  0.3,  q  =  80°  and  a  range  of  wave 
frequencies  from  w/fi  =  1.09  to  7.  A  typical  route  to  global  chaos  in  this  parameter  regime  is 
illustrated  in  Fig.  2.  For  c  =  0.14  which  is  below  the  stochasticity  threshold,  the  particle  orbits 
are  integrable  and  the  gain  in  energy  is  limited.  The  7m»x  oscillates  as  a  function  of  the  initial 
phase  (Fig.  2a)  with  different  phases  corresponding  to  different  values  of  H .  The  dependence  of 
the  on  the  initial  phase  means  that  the  stochasticity  threshold  is  also  dependent  on  the  phase. 
As  the  wave  amplitude  is  increased  to  c  =  0.155,  a  very  small  fraction  of  particles  are  stochsistically 
accelerated  (Fig.  2b).  At  e  =  0.27,  more  than  90%  of  the  particles  are  stochastic  (Fig.  2d).  The 
remaining  10%  of  the  particles  are  locally  stochastic  (weakly  chaotic)  but  are  limited  by  KAM 
surfaces  to  low  energies.  There  are,  however,  a  small  fraction  of  particles  ~  5%  near  th  =  t/2  that 
remain  integrable  even  after  the  wave  amplitude  has  been  increased  to  c  =  0.5.  The  7mix  is  in  fact 
smallest  a,t  rp  =  ir/2  as  evident  in  Fig.  2a. 

It  is  important  to  note  that  particles  can  gain  energy  without  being  trapped.  The  energy  gain 
of  a  particle  starting  with  zero  energy  is  determined  by  how  far  the  particle  is  from  a  resonance 
and  how  strong  that  resonance  is.  The  different  initial  phases  correspond  to  different  H  values,  and 
the  distance  of  the  particle  from  the  resonance  as  well  as  the  strength  of  that  resonance,  changes 
as  a  function  of  H .  For  the  R-x  mode  and  the  parameters  of  this  section,  the  spread  in  H  due 
to  different  phases  is  small.  The  phase  effects  axe  much  more  important  for  L-o  mode  as  we  show 
shortly.  For  particles  pear  «>  =  tt/2,  the  particles  are  far  from  resonance  and  the  KAM  surfaces 
separating  the  particles  from  resonance  persist  to  high  wave  amplitudes. 

As  the  simulated  value  for  the  stochasticity  threshold  at  a  given  frequency,  we  choose  two 
amplitudes  ei%  and  €90%  which  correspond  to  amplitudes  for  which  1%  and  90%  of  the  particles 
first  become  stochastic,  respectively.  We  use  690%  rather  than  €ioo%  because  as  shown  in  Fig.  2. 
there  can  be  a  large  gap  between  £90%  and  fioo%j  while  the  gap  between  fi%  and  £90%  is  much 
smaller.  In  other  words,  only  90%  of  the  particles  are  easily  accelerated  in  the  parameter  regime 
above. 

.-\t  all  frequencies  in  Fig.  3.  except  u/'il  =  1.09,  the  particles  reach  energies  above  7  >  .3  in 
Q.t  =  1000  after  they  become  stochastic.  We  may  thus  determine  what  fraction  of  the  particles 
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are  stochcLStic.  Another  check  is  to  solve  the  orbits  up  to  ilt  -  3000  and  check  v'hether  they 
were  accelerated  beyond  7max  at  =  1000.  These  two  tests  yielded  similar  results  in  all  cases. 
For  ijj/'Q  =  1.09,  the  inde.x  of  refraction  is  small  and  even  the  stochastic  particles  could  not  be 
accelerated  easily  above  =  3.  At  this  frequency,  we  determined  that  7  >  2  is  a  good  condition 
for  stochasticity. 

The  comparison  of  the  analytical  and  simulated  values  of  stochasticity  threshold  is  shown  in 
Figs.  3a  and  3b.  The  theoretical  threshold  is  indicated  by  the  solid  curve.  The  stochasticity 
threshold  Cihr,  is  seen  to  increase  with  frequency.  The  value  of  does  not  vary  much,  whereas 
the  value  of  £90%  increases  rapidly  as  a  function  of  frequency.  The  analytical  threshold  is  in  better 
agreement  with  (1%  than  £90%  for  most  of  the  frequency  range  below  u/Q  ~  5  with  the  situation 
reversing  at  higher  frequencies.  The  overlap  condition  for  stochasticity  is  a  crude  and  simple  way 
of  estimating  the  condition  for  stochasticity  and  it  is  not  surprising  that  the  simulated  threshold 
values  deviate  from  the  analytical  results.  The  fact  that  the  wave  amplitudes  are  rather  large 
(fthrs  ~  OT)  also  contributes  to  the  error.  An  additional  factor  of  tt  in  the  stochasticity  threshold 
which  is  often  invoked  to  take  into  account  the  effects  of  higher  order  resonances  (Chirikov,  1979) 
does  not  seem  appropriate  for  our  case,  and  it  would  in  fact  make  the  differences  between  the 
analytical  and  simulated  values  more  drastic. 

At  low  frequencies,  the  gap  between  £1%  and  £90%  is  small  and  thus  the  transition  to  global 
chaos  is  more  abrupt  as  is  evident  from  Fig.  3c.  The  gap  widens  at  higher  frequencies  and  for 
u;/f2  ^  4  it  becomes  difficult  to  accelerate  more  than  80%  of  the  particles. 

Next,  we  consider  the  variation  of  £thr,  as  a  function  of  propagation  angle  a  for  the  same 
parameters  as  in  Fig. 3  and  w/fl  =  2.  Both  the  analytical  and  simulated  £  are  seen  in  Fig.  4  to 
first  drop  rapidly  and  then  flatten  oui  as  a  function  of  a.  The  rapid  drop  of  fihrs  is  not  surprising 
since  the  system  is  integrable  at  q  =  0°.  The  analytical  threshold  in  this  case  overestimates  the 
threshold  condition  by  roughly  a  factor  of  2.  The  analytical  estimate  for  the  overlap  of  the  first 
three  resonances  but  not  extending  to  7  =  1  are  slightly  below  £17^.  Had  we  replaced  the  factor  of  4 
in  ( 15)  by  a  factor  of  2t,  the  agreement  between  the  analytical  and  simulated  values  of  f  thrs  would 
have  much  better  for  the  parameters  of  Fig.  4.  Even  though  there  exist  theoretical  arguments  for 
the  presence  of  this  factor,  it  does  not  wmrk  for  all  frequencies. 
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(ii)  R-x  mode,  Up  /n  =  2 

Next,  we  consider  the  plasma  parameters  corresponding  to  the  dayside  ionosphere  wv/n  =  2. 
and  Q  =  80°.  For  this  parameter  regime  only  80%  of  particles  are  easily  accelerated  and  there 
exists  a  rather  large  gap  between  6^0%  a-nd  f9o%-  As  a  result,  we  have  plotted  Cso^  rather  than  c-.q-;; 
in  Fig.  5a.  The  value  of  7max  is  very  small  near  ^  =  7r/2  and  3x/2  and  particles  with  initial  phases 
close  to  these  values  remain  integrable  even  for  very  large  wave  amplitudes.  The  analytical  threshold 
lies  mostly  between  (i%  and  e8o%-  The  threshold  amplitudes  tend  to  be  higher  than  those  in  case 
(i).  This  is  mainly  due  to  the  smaller  indices  of  refraction.  The  number  of  overlapping  resonances 
and  the  trapping  width  are  strong  functions  of  N.  For  fV  <  1,  the  trapping  width  increases  as  a 
function  of  harmonic  number  up  to  some  critical  value  lent  beyond  which  it  drops  ofi  exponentially. 
The  smaller  N ,  the  smaller  E^rit  and  thus  the  fewer  the  number  of  overlapping  resonances.  Since  the 
index  of  refraction  at  a  given  frequency  is  smaller  for  Wp/fi  =  2  compared  to  that  for  Up/Q.  =  0.3, 
fewer  resonances  overlap  and  the  diffusion  is  weaker  for  w/fi  =  2.  Thus,  it  is  easier  to  accelerate 
electrons  in  the  nightside  ionosphere. 

(iii)  L-o  mode,  uip/Cl  =  0.3 

Figure  6  illustrates  the  transition  to  chaos  for  w/fi  =  2.6  and  q  =  80°.  .-Vs  mentioned  earlier, 
the  phase  dependence  of  the  is  more  important  for  the  L-o  mode  than  it  is  for  the  R-x  mode. 
This  is  particularly  true  at  large  angles  where  the  term  €(simh)/iV||  in  the  Hamiltonian  becomes 
significant.  In  Fig.  6a  all  the  orbits  are  integrable  and  the  minimum  acceleration  is  seen  to  occur 
at  ti’  =  0,  TT  and  27r.  As  the  wave  amplitude  is  changed  from  0.35  to  0.5.  only  ~  60%  of  the  particles 
become  stochastic.  The  particles  between  V'  ~  0  to  tt  remain  mostly  integrable.  It  is  not  until  £ 
has  been  raised  to  ~  1  that  ~  50%  of  the  particles  between  th  ~  0  — '  tt  become  stochastic.  The 
phases  between  ~  0  to  t  correspond  to  positive  H  values.  Recalling  that  p:|  a  {tQ/u  -  H)  and 
A7  a  it  is  easily  seen  that  the  trapping  width  is  larger  and  the  Cihrs  smaller  for  negative  H 

vedues. 

The  behavior  shown  in  Fig.  6  is  typical  of  what  happens  at  other  frequencies  in  this  parameter 
regime.  Thus,  we  use  (1%  and  (gQ%  as  the  experimental  estimates  for  Cihrs-  The  comparison  of 
the  analytical  and  experimental  values  of  ft(,„  are  shown  in  Fig.  7a.  The  analytical  value  of 
changes  dramatically  as  a  function  of  frequency  whereas  the  experimental  value.s  vary  slowly  witii 
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frequency.  The  analytical  value  of  fthrs  is  seen  to  be  highly  inaccurate  for  <  2  and  ujil  >  5 
due  to  the  importance  of  the  phase  effects  for  the  L-o  mode.  If  we  replace  by  //  and  by  the 
canonical  momentum  in  (15).  we  find  much  better  agreement  with  the  experimental  values  as  we 
will  show  shortly. 

The  ratio  of  f6o%  to  ti%  fluctuates  between  ~  1.15  to  2.15  as  evident  from  Fig.  7b.  This  ratio 
gets  smaller  at  higher  frequencies  ir.  contrast  to  case  (f).  The  fthrs  values  are  typically  larger  than 
those  for  the  R-x  mode  in  Fig.  3  in  agreement  with  the  previous  predictions  [Karimabadi  et  al.. 
1990]. 

Finally,  the  variation  of  fthrs  with  the  propagation  angle  a  is  shown  in  Fig.  8.  The  variation 
of  fthrs  with  angle  has  the  general  form  of  that  for  the  R-x  mode.  The  theoretical  estimates  arc 
in  reasonable  agreement  with  experimental  results  for  a  between  10°  and  80°.  For  a  >  80°.  the 
theoretical  threshold  (15)  appears  to  be  inaccurate.  As  we  mentioned  e’^'ller,  this  is  due  to  the  fact 
that  at  large  angles  the  phase  dependent  term  in  R  due  to  L-o  mode  becomes  very  large.  There  is 
a  simple  way  of  modifying  (15)  to  account  for  the  phase  dependence:  replace  Ho  by  H  and  by 
the  canonical  momentum.  Thus,  depending  on  the  initial  phase,  we  would  have  a  different  v.nlue  of 
//  which  we  then  use  in  (15).  The  accuracy  of  this  procedure  is  tested  in  Fig.  9.  We  find  a  good 
agreement  with  the  experimental  result.  The  agreement  is  surprisingly  good  considering  the  crude 
way  that  we  have  incorporated  the  phase  effects  into  the  expression  for  the  trapping  width. 

The  acceleration  of  particles  starting  at  rest  generally  requires  larger  '  ave  amplitudes  than 
those  starting  with  a  finite  energy.  The  Bessel  functions  in  the  trapping  width  go  to  zero  if  evaluated 
at  7  =  1.  and  there  are  usually  persistent  KAM  surfaces  that  separate  such  particles  from  the  nearby 
resonance.  The  analytical  results  are  however  in  reasonably  good  agreement  with  the  experimental 
results  in  spite  of  the  complexity  of  the  dynamics  at  low  energies. 

VT  Test  of  the  Quasilinear  Theory^ 

In  order  to  test  the  quasilinear  theory,  we  have  solved  the  diffusion  equation  for  many  different 
parameters  and  compared  the  results  with  e.vact  orbit  calculations.  We  take  the  initial  particle 
distribution  function  to  be  a  delta  function  centered  at  7=1.  In  the  exact  orbit  calculations,  we 
follow  the  orbits  of  256  particles.  <ali  with  initial  7  =  1  and  phases  distributed  uniformly  between 
0  and  2t.  up  tn  []t  -  -iOOO.  This  corresponds  rou.shly  to  the  relevant  length  scale  for  ionospheric 
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heating  [Menyuk.  et  al..  1988].  The  resulting  distribution  function  is  very  bumpy  due  to  the  sn.dii 
number  of  particles  used.  To  obtain  smooth  distributions,  many  more  particles  wouid  have  to 
be  used — an  approach  which  is  both  both  impractical  and  unnicessary.  Instead,  we  compare  the 
moments  of  the  two  distribution  functions.  The  small  number  of  particles  which  we  use  is  suffic.e.nt 
in  this  case  to  yield  an  accurate  result.  In  the  following,  we  use  the  quantities  (7  -  7o)  =  (Al) 

{(7  ~  lo]')  =  {(^7)')  to  make  the  comparison,  where  70  is  the  initial  7  of  the  particles  and  the 
brackets  indicate  an  average  over  the  initial  conditions: 


((7  -  lo)  )  = 


jrii-'io)^f{i,t)di 


/r/(7.fy7 

Note  that  D-y-y  a  ((A7)*).  The  tes’’  of  the  quasilinear  theory  is  presented  below. 
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(i)  R-x  mode,  u/Q,  =  1.97 

First,  we  consider  the  diffusion  due  to  a  R-x  mode  with  u>/Q  =  1.97.  q  =  80°  and  Wp/fl  =  0.3. 
The  value  u/p/n  =  0.3  corresponds  roughly  to  the  nighttime  ionosphere  at  130  km  [Gurevich,  1978]. 
The  theoretical  and  numerical  Cihr*  are  0.19  and  0.143  respectively. 

We  carry  out  our  comparison  for  wave  amplitudes  ranging  from  0.19  to  10.  The  variation  of 
several  quantities  as  functions  of  time  are  shown  in  Fig.  10  for  c  =  G.19.  The  quantity  ((A7)") 
is  seen  to  follow  a  linear  increase  in  time  in  accord  with  the  underlying  assumptions  built  into 
the  Fokker- Planck  equation.  The  particles  have  not  had  enough  time  to  sample  all  the  available 
resonances  and  as  a  result  the  quantity  ((A7)*)  increases  with  time.  Once  the  particles  reach  the 
energies  beyond  which  D-yy  becomes  very  small,  the  ((Ay)')  flattens  out  as  a  function  of  time  and 
diffusion  stops.  For  the  parameters  of  Fig.  10,  the  flattening  out  does  not  happen  until  Qt  >  3000. 

The  quantity  DIDqi  which  is  the  ratio  of  the  exact  orbit  calculation  of  ((A7)'i  to  that 
obtained  from  the  diffusion  code,  deviates  significantly  from  unity  at  early  times  (Fig.  10c'  but 
settles  down  close  to  its  asymptotic  value  for  ill  ^  750.  This  is  expected  since  from  Fnis.  lf)b  and 
IGc.  it  follows  tbn,  the  diffusion  formalism  is  strictly  valid  for  times  long  corapareii  with  corroiaf.on 
time  scales.  Simihar  results  arc  also  seen  at  higher  amplitudes  ’is  is  apparent  in  i  lgs.  !1  and  id. 

The  deviation  of  particle  .iiotion  from  the  predicted  quasiiinear  diffus.on  at  eriri\'  times  cmi 
oc^'iir  due  to  at  least  t'vo  etTevis:  First,  there  can  exist  smail  but  finite  ^iand.s  01  stability 
an  otherwise  stocha.stic  phase  space.  .-V  pnrti.-le  <'  .ning  close  to  such  an  :,d;and  ran  get  tempore-  :. 
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trapped.  Such  stickiness  in  phase  space  can  obviously  have  an  effect  on  the  diffusion  of  the  pariicies. 
As  we  mentioned  earlier,  tlie  trapping  width  increases  as  a  function  of  harmonic  number  ana  in 
the  transition  to  global  chao.s,  the  border  of  chaos  lies  at  the  lowest  energies.  Since  w’e  start  the 
particles  at  7  =  1.  and  the  stickiness  is  most  important  at  the  border  of  chaos  which  is  also  at 
7  ~  i.  the  particles  often  experience  some  initial  stickiness;  however,  once  the  particles  reach  the 
higher  energies,  they  diffuse  freely.  Thus,  there  exists  a  finite  time  before  all  the  particles  can  diffu.se 
freely.  The  second  source  of  deviation  from  the  predicted  quasilinear  behavior  is  that  the  particles 
may  initially  sample  a  few  resonances  in  which  case  the  motion  is  closer  to  a  coherent  acceleration 
than  diffusion.  The  first  effect  results  in  a  retardation  in  the  diffusion  and  the  second  effect  leads 
to  diffuaion  larger  than  Dqr^.  The  balance  between  these  two  effects  determines  D/Dqi  at  early 
limes.  For  il.t  <;  7.50,  the  above  two  effects  become  much  less  important. 

The  Fig.  lOd-e,  show  the  evolution  of  {z)  and  (z)_^  >  in  time,  where  {:)  is  the  average 
distance  in  the  2  — direction  iravelied  by  the  particles  whereas  {z)_^  >  is  the  average  z-distance 
travelled  by  only  those  particles  that  attained  7max  ~  10  during  the  length  of  the  run.  As  the  wave 
amplitude  is  increased,  more  of  the  particles  can  reach  p'max  =  10  and  the  evolution  of  (z)  and 
(z)_^  >  10  I'^come  more  similar  as  shown  in  Figs.  1’  and  12.  For  parameters  of  this  section,  nearly 
all  of  the  particles  arc  accelerated  in  the  positi^•e  z-direction.  We  can  estimate  the  variation  of  z 
with  time  from  (2);  7  -  yc./fcA’ii)  ~  1  or 

r.  =  i"  -  li,Vi|c/-.  !29) 
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The  fraction  of  particles  having  •>  >  10  is  plotted  in  Fig.  lOe.  The  theoretical  curve  foilov.s 
the  e.xperimental  curv<;'  very  closely  especially  at  later  times.  The  same  i.s  al.so  true  at  oilier  ■>'  .ave 
amplitudes  as  shown  in  Figs.  11  and  12. 

N'e.xt  we  plot  values  of  {A-f) / (^i) q and  Df  Dqc  at  fit  =  3000  for  a  ranee  of  wave  anipiitudes 
from  f  =  0.19  to  10.  Here,  (A7}  and  (1^7) ql  are  the  average  (over  the  initial  condition  1  ol  A- 
as  obtained  from  numerical  solutions  of  orbits  and  the  diffusion  code,  respectively.  The  results 
are  shown  in  Fig.  13.  The  quasiiinear  theory  appears  to  be  quite  adequate  in  describing  the 
pajticle  motion  over  this  rather  large  range  in  wave  amplitude.  For  c  10.  the  particle  motion 
approaches  the  integrable  regime  of  unmagnetized  plasma  in  which  case  the  diffusion  equation 
would  no  longer  be  valid.  The  deviation  of  D/Dqi  from  unity  as  a  funcliun  of  wave  amplitude 
has  a  sporadic  behavior  and  does  not  seem  to  follow  a  nice  oscillating  pattern.  In  short,  wc  have 
found  no  evidence  for  oscillations  of  D/Dql  a-  function  of  wave  amplitude,  in  contrast  to  what 
was  found  by  Rechester  and  White  [1980]  in  their  studies  of  the  standard  mapping.  Furthermore, 
the  fluctuation  amplitude  of  DjDqi  is  much  smaller  than  that  seen  in  their  work.  The  presence  of 
resonances  of  different  sizes  clearly  leads  to  averaging  beyond  that  which  is  found  in  the  standard 
map  and  thus  to  a  closer  correspondence  to  quasiiinear  theory.  These  results  suggest  that  the 
standard  map  is  not  am  ideal  paradigm  for  real,  physical  systems. 

(ii)  L-o  mode,  ui/fl  =  2.6 

Next,  we  consider  the  diffusion  due  to  the  L-o  mode  with  w/Tl  =  2.6,  q  =  SO®  and  u-p/ll  =  i).3. 

The  theoretical  threshold  is  0.42  and  the  numerical  thresholds  fi%  and  (^0%  O  O.f.T 

lespectively. 

The  evolution  of  the  various  quantities  as  functions  of  time  are  shown  in  I  igb.  i  1  throueii 
The  quantity  ((iNy)^)  follows  a  more  or  less  linear  increase  in  time  as  before. 

The  diffusion  is  larger  than  Dqi  at  early  times  as  shown  in  Fig.  14c  even  though  for  f  =  ;27 

U’ss  than  .60%  of  the  particles  are  stochastic,  bimilar  behavior  is  also  seen  at  higher  v.uue  ampiiluoes 
(Figs.  Me.  16c).  This  is  mainly  due  to  the  fact  that  at  early  times  some  particles  a.'c  coherenii\’ 
acceierated  by  a  few  resonances.  .As  before,  for  ilt  >  7.30.  the  D{Dq!,  .settles  cio'.vn  close  tu  i:.- 

final  value.  For  (  —  0.427.  the  asymptotic  value  of  DIDql  is  below  unity  (Fig.  1  Ici.  1  hi-s  n  ii'ie 

to  the  fact  that  at  this  wave  amplitude  less  than  60%  of  the  particles  are  stochastic.  .As  the  w.ive 
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amplitude  is  increased,  more  of  the  particles  become  stochastic  and  eventually  D/ Dql  becomes 
larger  than  unity  as  shown  in  Figs.  15c  and  IGc. 


The  evolution  of  {z)  and  (’)_  >  in  time  are  shown  in  Fig.  Md  for  e  =  0.427.  .-V  good 
majority  of  the  particles  are  now  travelling  in  the  negative  rr-direction  and  the  analytical  estimate 
(30)  agrees  reasonably  well  with  but  not  with  (7).  T"his  is  expected  since  for  c  =  0.427 

less  than  60%  of  the  particles  are  stochastic  and  f30)  is  strictlv  valid  for  high  energy  particles. 
Furthermore,  as  we  have  already  shown,  for  the  L-o  inode  the  phase  effects  become  important  and 
in  estimating  t’..,  it  is  the  full  Hamiltonian  //  and  not  //„  that  has  to  be  used: 

A'li 


V,  =  (y  —  /{ jmc  -f  — -rm  sin  a  sin  tg) 

l'?l  ^!! 


131) 


N'ow.  V,  can  be  either  negative  or  positive  depending  on  v  and  II .  Notice  that  (z)  >  is  actually 

larger  in  Fig.  i4e  than  it  is  in  Figs.  15o  and  IGe.  This  is  a  result  of  the  high  energy  particles 
moving  in  the  positive  ^-direction  when  t  =  0.427:  whereas,  when  f  =  3.9  and  10,  the  high  energy 
particles  move  large  distances  in  both  the  positive  and  negative  z— directions,  and.  thus,  averaging 
yields  smaller  values  j-vaiues.  When  e  =  0.427,  the  particles  having  7^10  cover  a  range  in  ;  from 
0  to  27  km  in  Qt  =  3000;  whereas,  for  c  —  3-9  and  10.  the  range  in  z  is  —25  to  40  km  and  -30  to 
43  kra.  respectively. 

The  fraction  of  particles  with  7  >  10  as  obtained  from  the  diffusion  code  and  exact  orbit 
calculations  are  shown  in  Figs.  14e,  I'e  and  IGe.  The  agreement  between  the  two  methods  seems 
very  good  for  Qt  >  1000. 

Finally,  we  examine  the  behavior  of  i Ay) / {S.')  q i  and  D/ Dqi  as  functions  of  wa.ve  amplitude. 
This  is  shown  in  Fig. 17  where  we  have  nioitoa  {A',)  j  [Ay)  qi,  and  D/ Dqi,  at  ilt  =  3000  for  a  ranee 
of  wave  amplitudes.  I  he  quantity  (.A7)/ (Ay)  is  well  below  unity  for  t  =  0.43  and  it  docs  :ioi. 
reach  unity  until  the  wave  .amplitude  h:i.s  been  increased  to  c  ~  1.5.  For  t  21  1.5,  ^hc  (Ay)/{A-  )  qi 
has  small  fluctuations  aiiout  nnitv.  The  {Ay)Q[_  was  obtained  by  using  H  =  i.  However,  a.s  we 
-iioweri  earlier,  the  iilia.sc  etu-cts  ar->  important  for  the  l,.<j  mode  and  the  tihrs  depends  stroneiv  or, 
the  vaiue  of  II  which  h.a.s  a  unite  range  ’v.-n  for  7  1;  I .  The  t.sr,  can  change  by  more  than  .a  facto.r 
of  3  depending  on  the  v.alue  of  II  (see  Fi.T.  qi.  'Fl-us.  wc  t.rke  into  arcrunt  tiic  fraction 

o!  particles  that  are  not  sto(.ha.--tic  at  c  given  wave  amplitude.  The  ratio  of  .stocha.'tic  particles  i,, 
-  G0%,  at  ,  3:. 


ii. ,'5%  at  (  —  O.'h.  "Oyc  at  '  —  i.l  03%  .it,  c  -  1.3.  and  '19%  at  c  -  1  5. 


This  is  why  (A~i) / {Ay) is  below  unity  for  e  ^  i.5.  Let  us  denote  the  integrable  particles  by  tix- 
subscript  ‘iiit’  zind  the  stochastic  particles  by  the  subscript  ‘st’.  What  we  should  really  be  plotting  in 
Fig.  iTais  {Ay)  ,t / {Ay) q i  and  not  {Ay)/{Ay)QL.  But  {y),tl{l)QL  ~  (7)(-'V.nt  --  /gL- 

where  we  have  assumed  {y),ntl {1')ql  ^  1-  Using  the  values  of  N,t  cited  above  for  various  e's.  it 
then  follows  that  {Ay) ,t/ {Ay)QL  is  indeed  very  close  to  unity.  In  short,  the  discrepancy  between 
{Ay)  and  {Ay)QL  in  Fig.  17a  for  e  <  1.5,  results  from  the  inclusion  of  integrable  orbits  in  {y  )  but 
not  in  {Ay)Qi.  The  fraction  of  stochastic  particles  can  be  calculated  analytically  by  incorporating 
phase  effects  in  (15)  as  we  did  in  Fig.  9.  We  then  find  that  stochastic  particles  are  well  described 
by  the  diffusion  code. 

Similarly,  the  agreement  between  D  and  Dql  is  very  good.  This  result  is  encouraging,  consid¬ 
ering  that  the  theoretical  analysis  presented  in  this  paper  are  at  their  w'orst  for  the  parameters  of 
this  section  where  the  phase  effects  are  important  and  H  is  not  just  confined  to  one  value.  Thus, 
we  conclude  that  the  diffusion  formalism  is  quite  robust  and  is  highly  accurate  in  predicting  the 
time  evolution  of  an  ensemble  of  particles. 


VI.  Application  to  the  Ionosphere 

In  this  section,  we  tissess  the  plausibility  of  accelerating  the  electrons  in  the  ionosphere  using 
ground  based  transmitters. 

There  are  three  requirements  for  the  effective  acceleration  of  particles:  (1)  the  wave  amplitude 
must  be  larger  than  Cihrsi  (2)  particles  must  remain  in  the  system  f ionosphere)  long  enougii  to 
diffuse  to  high  energies;  and  (3)  the  acceleration  mechanism  should  be  insensitive  to  details  of  the 
initial  particle  distribution. 

If  condition  (1)  can  be  met,  then  condition  (3)  is  also  satisfied  since  the  stochastic  part:cie 
.acceleration,  in  contrast  to  any  coherent  mechanism,  is  quite  robust.  First,  we  consider  the  power 
needed  in  order  to  achieve  f  ^  tihrs  in  the  ionosphere.  As  we  showed  earlier,  the  lower'  o  r  irs 

for  the  R-x  mode  in  the  nighttime  ionosphere  (wp/fl  =  0.3)  in  the  frequency  range  19  -  2. 

The  threshold  can  be  as  low  as  0.08  at  a  ~  50°  and  r  0.1  at  o  ~  80°.  At  such  low  wave 
amplitudes,  only  a  small  fraction  Z%  of  particles  can  be  accelerated.  Fortunately,  the  transition 
to  global  chaos  is  abrupt  in  this  parameter  regime  and  once  c  is  raised  tc  ~  0.14,  more  than  90% 
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t>f  the  particles  become  stochastic.  Thus,  the  ininiruum  wave  amplitude  for  stochastic  particle 
acceleration  by  waves  with  ;V||  <  1  is  between  t  ~  (i.OS  -  0.1-1. 

The  stochasticiiy  tlireshold  is  in  general  much  lower  (c  —  in'"*  -  10“^)  for  waves  having  .'Vi.  >  1, 
but  the  Hamiltoniart  surfaces  are  closed  in  such  cases,  only  a  few  resonances  are  available,  and  the 
acceleration  is  much  smaller  than  in  cases  where  !V.|  <  1. 

The  power  flux  is  related  to  the  wave  amplitude  by  [Menyuk,  et  ai,  19S8]: 


Thus,  a  power  flux  in  the  range  ~  0.7  —  2. .3  W/cm~  is  required  to  accelerate  the  zero  energy 
electrons  to  large  energies  in  the  ionosphere.  Let  us  for  the  moment  suppose  that  such  a  pow.er 
‘lux  can  be  achieved.  \\’e  then  e.stimate  the  size  of  the  acceleration  region.  The  main  motion  of 
the  particle  is  in  the  ^-direction  and  the  particles  e.xecute  Larmor  motion  in  the  perpendicular 
motion  with  r [^  —  v±~i/Q.  Using  (30)  and  setting  z  ~  ri,  we  obtain  7t’j_/c  ~  /ViifiT  For  the 
parameters  in  Fig.  10  and  for  =  3000,  we  have  7t’x/c  ~  7  ~  512.  Thus,  the  escape  of  the 
particles  occurs  in  the  z-direction.  Even  though  the  evolution  of  <  z  in  time  is  insensitive 

to  the  wav.e  amplitude,  more  particles  can  be  accelerated  in  less  time  at  larger  wave  amplitudes. 
For  e  'z;  0.1  —  0.2,  the  typical  distance  that  is  required  for  accelerating  a  significant  number  of 
particles  to  7  >  10  is  of  the  order  of  10  —  20  km  as  shown  in  Fig.  10.  Thus,  considering  a  region  of 
ionosphere  10  km  x  20  krn,  the  total  power  required  is  ~  1.4  x  10'"  -  4.6  x  10'"  W.  If  we  launcli  1 
msec  pulses  with  a  duty  cycle  of  1  per  minute,  we  find  an  average  power  of  ~  iO'  -  10'’  These 
values  are  beyond  the  present-day  capabilities. 

In  view'  of  cuir  lindinus,  we  are  forced  to  Cvjnsider  ways  to  reduce  the  stochasticity  threshold. 
There  are  several  possibilities;  fl)  Prcaccelcrate  the  electrons.  This  involves  a  two-step  process. 
First,  the  particle.s  are  accelerated  to  weakly  relativistic  energies  and  then  the  rt  waves  are  applied. 
Since  the  Cii,„  is  lower  at  higher  energies,  we  could  make  the  particles  stochastic  for  lower  rf 
amplitudes.  The  deperuience  of  Cthrj  on  energy  is,  iiowever,  weak  and  this  techiii(|ue  can  at  most 
lower  the  amplitude  required  by  a  factor  of  2.  (2)  .Apply  several  waves  at  closely  spaced  frequencies. 
The  problem  is  that  if  the  wave  amplitudes  are  too  ;  'l‘•i()w  tiie  co,rj  value  ot  a  single  wave.  i.ho 

resulting  diffusion  would  be  limited  and  slow.  (3)  1  he  magnetic  field  in  the  ionosphere  is  well 
represented  hv  magnetic  d'pole.  Thus,  a  more  realistic  application  to  the  ionosphere  requires 


taking  into  account  the  inhomogeneity  of  the  magnetic  field.  It  is  well  known  that  even  a  ■veak 
inhomogeneity  in  the  magnetic  field  can  lower  the  Cihrs  fiy  much  as  2  orders  of  magnitnce. 
Physically,  the  change  in  Bo  is  equivalent  to  having  several  resonances  for  a  given  harmonic  nr  ;n;,e:r 
Thus,  the  resonances  are  more  tightly  packed  and  the  overlap  occurs  at  smaller  amplitudes  timn 
in  the  uniform  field  limit.  In  order  to  have  a  strong  enough  diffusion,  however,  the  wave  ampiitude 
cannot  be  much  smaller  than  in  the  uniform  field  limit. 

Finally,  the  self-consistent  effects  are  expected  to  play  a  role  once  a  significant  number  of 
particles  are  accelerated.  The  study  of  self-consistent  effects  tis  well  as  mode  conversion  are  currently 
underway,  and  preliminary  results  were  reported  by  Akimoto  and  Karimabadi  (1989). 

VII.  Conclusions 

In  this  work,  we  have  carried  out  a  detailed  comparison  of  the  theoretical  stochasticity  threshold 
with  exact  orbit  calculations.  We  reduced  the  diffusion  equation  to  a  one-dimensional  form  by- 
transforming  to  a  coordinate  system  where  the  Hamiltonian  surface  is  one  of  the  axes.  We  then 
constructed  a  model  for  the  diffusion  coefficient  and  wrote  a  finite  element  code  that  solves  the 
diffusion  equation.  This  code  has  been  carefully  compared  with  a  particle  code,  and  it  shows 
excellent  agreement  over  a  wide  range  of  frequencies  and  wave  amplitudes.  We  applied  our  results 
to  the  problem  of  radio-frequency  acceleration  of  the  electrons  in  the  ionosphere.  We  found  that 
the  powers  needed  are  beyond  present-day  technology.  We  then  discussed  several  ways  to  overcome 
this  problem. 

Finally,  we  emphasize  that  the  method  of  reducing  the  diffusion  equation  to  one  dimension  and 
the  code  that  we  have  developed  are  quite  powerful  and  are  expected  to  find  many  applicat.ons. 
We  are  currently  applying  the  above  techniques  to  the  problem  of  particle  acceleration  in  planetary- 
shocks. 

.Acknowled,  ment 

The  authors  would  like  to  thank  Dr.  E.  Moghadam  Taaheri  for  useful  discussions  concerning 
the  finite-element  technique.  This  work  was  supported  by  the  Space  Physics  Theory  ITogram  at 
rC  .San  Diego  and  tlie  .Air  F<  rr  Geophysics  Laboratory  at  Hanscoiii  .Air  Force  Base.  Gompiitine 
was  p<M formed  on  the  San  Dn-eo  .7upe. computer  Center  CR.AA’  A'-MP. 


8k 


Appendix  A:  Diffusion  Equation  Along  the  H-surface 

In  this  appendix,  we  show  how  the  diffusion  equation  can  be  reduced  to  a  one-dimensionai 
form  by  choosing  a  coordinate  system  that  has  H  as  one  of  the  axes.  A  given  // -surface  is  given 
by; 

Ho  =  ’nc--)o  - 


W’e  can  rewrite  the  equation  for  Ho  as 


where 


Pj 

m^c- 


1  - 


(A  •2a) 


bn  (A2c) 

A II 

The  H -surfaces  are  not  confocal  since  both  an  and  6//  depend  on  H.  Thus,  a  confocal  orthogonal 
transformation  cannot  be  made.  In  fact  there  is  no  simple  transformation  from  (P||,Pi)  to  an 
orthogonal  coordinate  system  with  the  //-surface  as  one  of  the  axes.  Fortunately,  the  orthogonality 
condition  is  not  necessary,  and,  as  we  now  show,  any  component  of  the  diffusion  tensor  not  parallel 
to  the  //-surface  is  zero. 

Let  us  consider  the  diffusion  equation  under  the  transformation  from  (py.pi.)  to  (H,g)  where 
fj  is  arbitrary  for  the  moment.  We  start  with  (19)  rewritten  in  tenso.  notation: 


dl  _ 

dt  R  (Jp^ 


(A.3) 


where  /?  =  pi  is  the  .lacobian  of  the  transformation  from  Cartesian  to  cylindrical  coordiiiaios. 


From  the  chain  rule,  we  have 

d  _  dyk  d 
Op,  Op,  Ovk ' 

where  -  //  V.’  “  .d-  L'sing  this  in  (.A3),  yields 


(A4) 


‘II  l^JL  (  HD  I]iL^\ 

Ot  R  Op,  Opk  \  '-OpjOy,) 

0  /  Ovk  -  I]H.IIJ\  -?n  dyi  f  _2_  ( }_I]il^\ 

Oyk\dp,  '’OpjdytJ  ‘  ’’ Op;  dyt  dyk  \  R  dp,  )  ' 
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or 


^  =  IL]  ^  RD 

dt  Oijk  V  0y()  ''Op,  dye  Oy^  \Rdp,  J  ' 


where  Dke  is  the  new  diffusion  tensor  and  is  given  by 

^  _dyk  dye 
L>kt  —  -rt — ^•’'3 — 
dpi  '  dpj 


Writing  out  the  various  components  of  the  diffusion  tensor  explicitly,  we  have 

^  dHdH 


r  /  /)  //  \  ^ 

Dhh  =  ^[d^Uu[^^J  +2^x11, 


(ASa) 


^39  ~ 


(  d9\  „  ^  (  dg 

V^Pll  j  “  9p||  Spx  \5px 


i  ASb) 


Dng  =  5] 


5p||  dp|| 


-LX< 


gif  gg 

dpL  dpi.  J  ’ 


dg  dH  \ 

V5p||  opx  5p||  ^pxj  ''' 


and 


Dng  =  DgR. 


(A8c) 


(ASd) 


The  derivatives  are  evaluated  at  the  crossing  points  of  a  given  if —surface  and  the  resonanc"  curves, 
fn  particular,  we  find 

dH  n  I  dH 

- = - .  (  \f'  I 

5pii  w  pxA'ii  5px 

It  follows  that  Dhh  is  identically  zero.  Furthermore,  using  the  relation  and  substituting  '  22) 
in  (A8c),  we  find  that  D^g  is  identically  zero  independent  of  the  form  of  g.  So,  the  new  coordinate 
system  need  not  be  orthogonal.  The  fact  that  Dng  =  Dhh  =  0  is  hardly  surprising  since  both 
Dng  and  Z?////  involve  integrals  over  the  quantity  H  which  is  zero. 

A  useful  choice  for  g  is  the  relativistic  factor  7.  The  surfaces  of  constant  -  define  circles  in 
momentum  space.  The  o'd  variables  are  realted  to  the  new  variables  by 


me 


=  (1 - ^)A'l| 

me' 


(AiO.:, 


Px 

me 


7-  -  1  -  iV||^(7  - 


H 


me' 
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(A  10b  I 


.iiid  the  Jacobian  of  the  transformation  in  dimensionless  units  is  given  by 


R  —  "f  A II ,  (All ) 

which  is  valid  for  a  ^  90°.  Replacing  g  by  7,  and  using  (A8)  through  (All)  in  (A6)  and  (A7).  we 
obtain  the  diffusion  equation  (23a)  and  diffusion  coefficient  (23b). 
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